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ABSTRACT 


Traditional  methods  of  time  difference  of  arrival  (TDOA)  determination  suffer 
significantly  in  environments  fraught  with  co-channel  interference  and  low  signal  to  noise 
ratios.  Cyclostationary  signal  processing  techniques  offer  solutions  to  the  shortcomings 
of  the  traditional  TDOA  methods.  Specifically,  the  Spectral  Coherence  Alignment 
(SPECCOA)  method  of  TDOA  determination,  developed  by  the  Mission  Research  Corp. 
and  Statistical  Signal  Processing  Inc.,  performs  exceptionally  in  very  poor  signal  to  noise 
ratio  environments.  The  Applied  Research  Lab  at  the  University  of  Texas  at  Austin 
(ARL:UT)  has  developed  a  prototype  TDOA  system,  the  Carry-on  Multi-platform  GPS 
Assisted  Time  Difference  of  Arrival  System  for  the  Naval  Information  Warfare  Activity. 

It  currently  utilizes  a  traditional  complex  ambiguity  function  (CAP)  to  determine  the 
TDOA(s)  between  multiple  observers  and  an  ARL:UT  developed  closed  form  solution  for 
the  geolocation  of  the  emitter.  The  work  presented  here  takes  the  first  step  in  applying 
SPECCOA  to  the  ARL:UT  system.  Coding  both  SPECCOA  and  the  ARL:UT  closed 
form  solution  in  MATLAB®  makes  possible  a  quantitative  comparison  between  the  CAP 
and  SPECCOA  using  ARL:UT  real  world  test  signals. 
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I. 


INTRODUCTION 


A.  BACKGROUND 

The  passive  geolocation  of  electromagnetic  (EM)  emitters  plays  an  increasingly 
important  role  in  all  three  aspects  of  Information  Warfare  (IW).  In  order  to  protect 
friendly  communications  from  hostile  jamming  and  interference,  information  warriors 
must  first  locate  the  source  of  that  jamming  and  interference.  The  geolocation  of  an 
enemy  radar  exploits  the  enemy  by  gathering  information  ascertained  from  his  own 
information  gathering  system.  Finally,  in  order  to  attack  an  enemy  communications  node, 
its  coordinates  must  be  passed  to  the  strike  aircraft  or  the  cruise  missile.  In  addition  to  its 
importance  in  the  military  arena  of  IW,  passive  geolocation  of  EM  emitters  finds  use  in 
law  enforcement  surveillance,  search  and  rescue  operations,  navigation  and  the 
enforcement  of  communications  regulations. 

Three  aspects  of  EM  waves  lend  themselves  to  exploitation  for  geolocation 
purposes.  The  path  of  an  EM  wave  can  be  closely  approximately  from  transmitter  to 
receiver  given  the  frequency  of  the  transmission  and  some  basic  characteristics  about  the 
atmosphere  through  which  it  traveled.  In  addition,  when  EM  waves  arrive  at  two  moving, 
spatially  separated  receivers,  the  receivers  measure  different  Doppler  shifts  in  the 
frequency  of  the  transmission.  Finally,  the  time  that  an  EM  wave  arrived  at  two  spatially 
separate  receivers  yields  valuable  information.  These  three  characteristics  naturally 
enough,  are  the  basis  for  the  three  basic  methods  of  geolocation.  Angle  of  Arrival  (AOA), 
Frequency  Difference  of  Arrival  (FDOA)  and  Time  Difference  of  Arrival  (TDOA). 

AOA  geolocation  involves  the  use  of  highly  directional  antenna  arrays.  By 
measuring  the  phase  difference  in  the  EM  wave  at  each  antenna  element  the  receiver  can 
calculate  the  direction  from  which  the  wave  emanated.  Drawing  a  line  from  the  precise 
position  of  the  receiver  in  this  direction  yields  a  vector  containing  the  position  of  the 
transmitter  called  a  line  of  bearing  (LOB).  Moving  the  receiver  and  taking  another 
measurement  or  using  an  additional  receiver  located  elsewhere  produces  another  LOB. 

At  the  intersection  the  these  LOB(s)  is  the  transmitter.  AOA  is  the  oldest  of  the  three 
methods  and  has  been  known  by  many  names  throughout  the  years  of  its  use.  Beginning 
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just  prior  to  World  War  II,  the  name  most  commonly  associated  with  AOA  was 
triangulation.  With  small  measurement  errors,  the  three  vectors  had  a  finite  width  when 
drawn  and  thus  intersected  to  form  a  small  triangle  within  which  the  emitter  was  located. 
This  point  is  illustrated  in  Figure  1-1  below. 


Figure  1-1  -  Angle  of  Arrival  geolocation 

For  moving  collectors,  the  Doppler  effect  ensures  that  the  receiver  will  measure  a 
different  frequency  of  arrival  than  it  would  as  a  stationary  collector.  Given  two  collectors 
that  measure  two  different  frequencies  of  arrival  (FOA),  the  difference  of  these  FOAs 
(FDOA)  can  be  used  to  determine  the  position  of  the  emitter.  FDOA  produces  a  locus  of 
points  called  an  isodop  along  which  the  emitter  lies.  This  isodop  represents  all  the 
locations  from  which  the  EM  wave  could  emanate  and  produce  the  difference  in  Doppler 
shifts  measured  between  the  two  receivers.  Figure  1-2  shows  a  two  dimensional 
depiction  of  a  family  of  FDOA  contours  in  the  plane  of  the  paper.  Each  contour 
represents  a  specific  FDOA  between  the  observers.  The  emitter  could  lie  anywhere  along 
the  appropriate  contour,  and  thus  multiple  measurements  are  required  to  determine  its 
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location  precisely.  While  these  contours  can  be  approximated  as  two  dimensional  when 
the  emitter  is  on  the  surface  of  the  earth,  for  an  emitter  in  three  dimensions  they  would  be 


Figure  1-2  -  FDOA  contours 

surfaces.  A  receiver  wishing  to  use  Doppler  for  geolocation  purposes  must  possess  the 
ability  to  measure  frequency  more  precisely  than  the  smallest  Doppler  shift  expected. 
Otherwise,  the  small  Doppler  shifts  go  undetected  as  they  are  below  the  noise  in  the 
frequency  measurement.  This  has  been  the  limiting  factor  in  the  usage  of  FDOA  for 
precise  geolocation.  Measurement  of  frequency  precisely  is  much  more  difficult  and 
costly  than  the  measurement  of  time  precisely. 

The  precise  measurement  of  time  has  been  revolutionized  by  the  high  speed 
digital  technology  of  the  Global  Positioning  System  (GPS)  and  synchronization  circuits 
for  atomic  time  standards.  Using  the  difference  in  arrival  times  of  an  EM  wave  at  two 
spatially  separated  receivers,  a  locus  of  points  called  an  isochron  includes  all  possible 


locations  from  which  the  EM  wave  could  emanate  and  arrive  at  the  two  receivers  at  the 
two  different  times  measured.  In  two  dimensions,  like  the  situation  that  approximately 


Figure  1-3  -  Generic  TDOA  estimator 

exists  when  an  emitter  is  constrained  to  the  surface  of  the  earth,  the  isochron  can 
be  represented  by  a  hyperbola.  In  three  dimensions,  the  hyperbola  becomes  a  hyperboloid 

of  revolution.  A  simple  block  diagram  for  a  generic  TDOA  estimator  appears  in  Figure 
1-3. 

The  traditional  methods  used  for  TDOA  determination  include  the  general  cross 
correlation  (GCC)  (in  the  absence  of  Doppler  shift)  and  the  cross  ambiguity  function 
(CAF)  (both  TDOA  and  FDOA  simultaneously).  These  two  methods  perform  well  in 
benign  signal  environments.  That  is,  with  a  signal  of  interest  (SOI)  with  a  high  signal  to 
noise  ratio  (SNR),  in  the  absence  of  jammers  and  other  signals  not  of  interest  (SNOI)  in 
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the  frequency  band  around  the  signal  of  interest  (SOI),  the  GCC  and  CAF  methods  are 
superb  estimators  of  TDOA.  However,  in  a  peacetime  EM  spectrum  increasingly  packed 
with  signals,  and  a  battlespace  spectrum  potentially  flooded  with  SNOI(s),  hostile 
jamming  and  weak  SOI(s)  the  GCC  and  CAF  methods  suffer  from  severely  decreased 
signal  to  noise  and  interference  ratio  (SNIR).  In  light  of  these  weaknesses  in  hostile 
environments,  a  relatively  new  concept  in  signal  processing  called  cyclostationarity  has 
been  applied  to  the  TDOA  problem. 

Cyclostationary  processes  are  characterized  by  time-periodic  statistics.  A  wide- 
sense  cyclostationary  signal  specifically  has  a  mean  and  an  autocorrelation  that  vary 
periodically  in  time.  These^periodicities  are  usually  hidden  in  traditional  second  order 
stationary  processing  techniques  such  as  the  autocorrelation  function  or  its  frequency 
domain  counterpart,  the  spectral  density  function.  However,  when  second  order 
cyclostationary  techniques  such  as  the  cyclic  autocorrelation  function  and  the  cyclic 
spectral  density  function  are  used,  wide  sense  cyclostationary  signals  readily  display  these 
hidden  periodicities.  These  periodicities  now  revealed  may  be  utilized  for  the  detection, 
classification  and  parameter  estimation  of  the  signal.  Among  those  parameters  that  may 
be  estimated  is  the  TDOA  of  a  signal  impinging  upon  two  spatially  separated  antennas. 

The  robustness  of  the  cyclostationary  technique  for  TDOA  lies  within  the  fact  that 
every  signal  has  a  unique  set  of  periodicities  that  depend  on  such  parameters  as 
modulation  scheme,  bit  rates  and  spreading  codes.  Two  signals  that  may  appear  as 
spectrally  overlapping  in  the  traditional  stationary  spectral  density  function  can  be 
separated  with  great  accuracy  by  the  cyclic  autocorrelation  function.  In  highly  corrupt 
environments,  cyclostationary  processing  techniques  provide  the  capability  to  select 
specific  signals  for  geolocation,  nearly  independent  of  the  presence  of  jamming,  inference 
and  SNOI(s)  that  may  spectrally  overlap  the  SOI. 

The  United  States  Navy  has  historically  relied  upon  AOA  techniques  for  direction 
finding.  Beginning  in  the  1950’s,  the  Navy  began  constructing  the  High  Frequency 
Direction  Finding  (HFDF)  network  of  circularly  disposed  antenna  arrays  (CDAA).  These 
CDAA(s)  consist  of  a  circular  array  of  elements  and  reflectors  that  serve  to  detect  the 
direction  of  incoming  energy.  Given  the  mass  of  recent  base  closures,  their  numbers  have 


5 


dwindled  significantly.  In  addition,  they  are  not  an  organic  asset  for  battle  group 
commanders  considering  the  current  tactical  communications  environment  that  relies  less 
on  long-haul  HF  and  more  on  Very  High  Frequency  (VHF)  and  Ultra  High 
Frequency(UHF)  point-to-point  links  and  Super  High  Frequency  (SHF)  satellite 
communications  links. 

In  addition  to  the  HFDF  network,  the  Navy  added  the  OUTBOARD  and 
OUTBOARD  II  systems  to  some  of  its  Spruance-class  destroyers.  These  too  are 
primarily  HFDF  assets  and  have  been  followed  by  the  current  construction  of  COMBAT 
DF  equipped  Essex-class  amphibious  ships  and  Arleigh  Burke-class  guided  missile 
destroyers.  These  three  shipboard  systems  have  proven  extremely  useful  during  the  past 
years  but  as  they  have  become  the  object  of  intense  scrutiny  with  shrinking  ship-building 
budgets  and  increasing  decommissionings.  Consequently,  the  Naval  Security  Group 
Support  Activity  (NSGSA),  in  June  1993,  contracted  with  the  Applied  Research  Lab  at 
the  University  of  Texas  at  Austin  (ARL:UT)  to  develop  an  affordable,  low-risk 
TDOA/FDOA  geolocation  system  using  commercial  off  the  shelf  (COTS)  and 
government  off  the  shelf  (GOTS)  technologies. 

The  Carry-on  Multi-platform  Global  Positioning  System  (GPS)  Assisted  TDOA 
System  was  tested  for  the  first  time  in  August  1995  after  fifteen  months  of  research  and 
development.  In  further  testing  off  the  coast  of  San  Diego,  California,  the  prototype 
system  geolocated  HF,  VHF  and  UHF  emitters  with  a  mean  square  error  of  approximately 
100  meters,  using  a  CAF  based  TDOA  determination  algorithm. 

As  noted  above,  the  CAF  has  performed  exceptionally  well  in  the  relatively 
favorable  conditions  during  prototype  testing  but  in  theory,  is  susceptible  to  low  SNR 
conditions  and  co-channel  interference.  Cyclostationary  processing  techniques  for  TDOA 
determination  offer  potential  improvement  to  the  performance  of  the  ARL:UT  system  in 
highly  corrupt  environments.  The  evaluation  of  a  cyclostationary  TDOA  algorithm  in 
conjunction  with  the  closed  form  geolocation  algorithm  used  in  the  ARL:UT  prototype  is 
the  primary  objective  of  the  work  presented  here. 

As  a  secondary  objective,  the  cyclostationary  TDOA  determination  algorithm  and 
the  closed  form  geolocation  algorithm  will  also  serve  as  the  first  two  blocks  for  testing  in 
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a  developing  geolocation  software  workbench  at  Naval  Postgraduate  School.  The  intent 
is  to  model  the  system  in  the  Simulink®  environment  developed  by  The  Mathworks  Inc. 
and  allow  different  algorithms  to  be  exchanged  by  merely  interchanging  the  appropriate 
Simulink®  blocks.  This  will  allow  for  the  testing  of  all  aspects  of  the  geolocation 
problem  from  data  conversion,  filtering,  AO  A,  FDOA  and  TDOA  determination  and 
various  geolocation  solutions. 

B.  ORGANIZATION 

The  traditional  TDOA  determination  algorithms,  the  GCC  and  the  CAP  along 
with  their  disadvantages  in  corrupt  environments  are  present  in  Chapter  II.  Following  the 
illustration  of  the  vulnerabilities  of  these  two  algorithms,  the  theory  of  cyclostationarity 
and  TDOA  determination  by  cyclostationary  techniques  are  presented  in  Chapter  EL 
Chapter  IV  depicts  the  closed  form  solution  to  the  TDOA  geolocation  problem  developed 
by  Dr.  Petre  Rusu  at  ARL:UT  for  the  prototype  geolocation  system  which  is  described  in 
Chapter  V. 

Chapters  VI,  VII,  VIE  and  IX  constitute  the  bulk  of  the  investigation, 
summarizing  the  new  algorithms  developed  in  MATLAB®,  the  test  plan  and  the  results 
of  those  tests  including  the  plans  for  Simulink®  block  development.  Finally,  Chapter  X 
draws  some  conclusions,  offers  explanations  for  the  anomalies  encountered  and  suggests 
areas  for  further  research. 
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II.  TRADITIONAL  TDOA  DETERMINATION 


A.  GENERALIZED  MODEL 

For  the  development  of  the  two  traditional  methods  of  TDOA  determination, 
assume  a  stationary  signal  from  a  remote  emitter  impinges  upon  two  spatially  separated 
antenna  elements  as  illustrated  in  Figure  2- 1 .  In  the  presence  of  additive  white  Gaussian 
noise  (AWGN)  noise,  the  two  signals  at  receivers  1  and  2  can  be  modeled  as 

x(t)  =  s(t)  +  nj(t)  (2-1) 

and 

y(t)  =  A-s(t-D)  +  n2(t)  (2-2) 

where  ni(t)  and  n2(t)  are  assumed  to  contain  only  AWGN  with  no  strong  SNOI(s)  in  the 
frequency  band  of  interest  for  the  SOI,  s(t).  A  represents  the  complex  valued  relative 
magnitude  and  phase  mismatch  between  the  receivers;  D  is  the  TDOA  of  the  SOI 
between  the  signals. 

It  follows  that  the  autocorrelation  and  cross-correlation  functions  are  given  by 


/?,(T)  =  i?,(T)  +  /?„,(T)  (2-3) 

i?/T)=|Ap/?,(T)  +  /?„^(T)  (2-4) 

R^A'^)  =  A-R,(x-D)  +  R„^^(x)  (2-5) 

and  the  respective  spectral  density  functions  are 

5,(/)  =  5,(/)  +  5„,(/)  (2-6) 

Syif)  =  \Afs^(f)  +  S„^(f)  (2-7) 
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+  (/) 


(2-8) 


These  relationships  contain  the  parameter  of  interest,  specifically,  D,  the  TDOA. 
The  next  two  sections  show  two  different  methods  used  to  extract  D  from  these 
equations.  The  first,  the  General  Cross  Correlation  (GCC)  function,  can  be  used  in  the 
absence  of  relative  motion  between  the  two  receivers  and  the  emitter.  The  second,  the 
Complex  Ambiguity  Function  (CAF),  can  be  used  to  estimate  FDOA  and  TDOA 
simultaneously  as  is  necessary  when  relative  motion  between  the  receivers  and  emitter 
exists. 


Figure  2-1  -  Signal  model  scenario 

Of  key  importance  to  the  development  illustrated  here  is  the  statistical 
independence  of  s{t),  n\{f)  and  niit)  and  the  absence  of  in  band  interference.  While  it  may 
be  argued  that  the  signals  and  noise  measured  at  the  two  receivers  can  be  correlated  to 
some  extent,  that  scenario  greatly  complicates  this  case  and  is  beyond  the  scope  of  this 
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demonstration.  The  derivations  that  follow  are  drawn  from  [1],  [2],  which  treat  the  GCC 
more  specifically  and  [3]  and  [4]  which  address  the  CAP  more  so  than  the  GCC.  All 
assume  the  three  components  to  be  uncorrelated. 

B.  GENERALIZED  CROSS  CORRELATION 

Applying  (2-5)  in  the  absence  of  SNOI  in  the  SOI  band  of  interest  and  statistically 
independent  s(f),  riiit)  and  n2{t)  yields 


/?,„(T)  =  A-i?,,(T-D)  +  /?„.„^(T) 

(2-9) 

This  function  will  peak  at  T=  D,  the  TDOA  between  receivers  1  and  2.  Because  ni{t)  and 

n2{t)  contain  merely  AWGN,  their  cross-correlation  term  in  (2-9)  above  reduces  the  SNR 

of  the  measurement  but  does  not  add  interference  in  the  form  of  SNOI(s).  Further 

analysis  using  (2-6)  -  (2-8)  in  a  similar  manner  reveals 

(2-10) 

S^if)  =  \Afs^if)  +  S^^(f) 

(2-11) 

(2-12) 

To  extract  D  from  (2-12),  first  note 

•s,(/)=r  -^'’^2 

=  0  otherwise 

(2-13) 

which  assumes  the  SOI  effectively  occupies  a  finite  bandwidth  B  around  the  carrier  or 
intermediate  frequency /o.  Taking  the  ratio  of  the  spectra  and  using  (2-10)  with  (2-13), 
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(2-14) 


SAf)  0 


/  -  —  </ +— 
Jo  2  *^  —Jo'  2 

otherwise 


and  taking  the  inverse  Fourier  transform  of  this  ratio  gives 


K(t:)=  J 

L-% 


SyAf) 

SAf) 


A  sin[7ufi('C-D)] 
Tt(x-D)l2 


cos[27t/„(i:-D)] 


which  peaks  at  T  =  Z).  This  in  turn  may  be  rewritten  as 


(2-15) 


h^(t)=  JW(/)5^,(/)e"^^'^#  (2-16) 


The  weighting  function,  W(f),  is  defined  in  this  instance  as  l/Sjc(f)  in  (2-16)  which 
is  the  best  choice  given  no  a  priori  information  about  the  SOI  [5].  In  addition,  it 
distinguishes  this  case  as  the  generalized  cross  correlation  method.  Assigning  a 
weighting  function,  W(f)  =  1  reduces  the  TDOA  determination  to  a  simple  cross¬ 
correlation  as  in  (2-5).  Given  prior  knowledge  of  the  noise  and  interference 
characteristics,  other  choices  for  W(f)  include  the  Roth  impulse,  SCOT  and  PHAT  among 
many  others  [3]  which  are  designed  to  reduce  specific  noise  and  interference  problems. 

It  is  clear  that,  in  the  absence  of  significant  noise,  co-channel  interference  and 
relative  motion  between  receivers  and  emitter,  the  GCC  produces  the  desired  estimate  of 
the  TDOA  of  a  signal  from  the  ratio  of  the  estimates  of  the  spectral  density  function  of 
the  signal  at  one  receiver  and  the  cross  spectral  density  function  of  the  two  received 
signals.  Figure  2-2  above  illustrates  the  generalized  cross  correlation  process.  The  two 
filters,  Hi(f)  and  H2(f)  are  specifically  designed  to  remove  out  of  band  interference. 
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While  in  practice  the  GCC  is  often  done  in  the  frequency  domain,  Figure  2-2  depicts  the 
time  domain  correlation  of  the  input  signals.  The  processes  are  theoretically  equivalent. 

Clearly,  the  presence  of  an  interference  signal  in  the  bandwidth  B  defining  the 
spectral  densities  would  corrupt  the  estimate.  In  addition,  as  is  shown  in  the  next  section, 
Doppler  shifts  at  each  receiver  also  must  be  accounted  for  in  order  to  accurately 
determine  TDOA. 


Estimate  of  TDOA 


Figure  2-2  -  GCC  block  diagram 


C.  COMPLEX  AMBIGUITY  FUNCTION 

The  complex  ambiguity  function  (CAF)  can  be  interpreted  as  an  extension  of  the 
GCC  for  moving  transmitters  and/or  receivers.  Stein  in  [4],  showed  that  in  order  to 
properly  determine  the  TDOA  between  two  receivers  in  the  presence  of  a  Doppler 
difference  at  each  receiver,  the  spectrum  of  one  of  the  signals  first  must  be  translated  in 
frequency  by  an  amount/equal  to  the  Doppler  difference  (FDOA)  measured  between  the 
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observers.  In  order  to  show  this,  a  Doppler  shift  is  introduced  into  the  generalized  model 
from  above.  Equations  (2-1)  and  (2-2)  can  now  be  written  as 

x(t)  =  s(t)+n,(t)  (2-17) 

yit)  =  As(t  -  -I-  (r)  (2- 1 8) 

where /d  is  the  Doppler  difference  as  measured  between  observer  1  and  2. 

Now,  in  place  of  calculating  TDOA  with  the  two  dimensional  GCC  from  (2-16), 
which  in  the  presence  of  significant  Doppler  could  peak  at  a  value  that  does  not  truly 
correspond  to  the  TDOA,  it  is  necessary  to  calculate  the  three  dimensional  CAP  given  by 

r 

T 

(2-19) 

0 

The  simultaneous  determination  of  the  TDOA  D  and  the  Doppler  shift /a  causes 
IA(D,/^)I  to  peak.  At_^  =  0,  the  CAP  reduces  to  a  GCC  problem  as  outlined  above.  Por 
fd^O,  the  CAP  may  be  thought  of  as  a  GCC  performed  after  frequency  shifting  the 
spectrum  of  y(t)  up  or  down  as  necessary  by  an  amount  equal  to /a.  A  block  diagram  of 
the  CAP  operation  appears  in  Figure  2-3.  Note  the  similarity  between  it  and  the  GCC 
diagram  of  Figure  2-1 . 

The  three-dimensional  width  of  the  correlation  lobe  is  directly  proportional  to  the 
accuracy  of  the  estimates  of  TDOA  and  FDOA.  Stein  further  points  out  in  [4]  that  the 
variance  of  the  estimate  for  each  parameter  can  be  related  to  the  noise  bandwidth,  the 
signal  bandwidth,  the  integration  time  and  the  effective  input  signal  to  noise  ratio  as 
follows 


(2-20) 
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Estimate  of  TDOA 


Figure  2-3  -  CAF  block  diagram 


_J _ 1_ 

~T,jBTy 


(2 


where 


B  =  noise  bandwidth  at  the  input  of  both  receivers 


P=27t 


|K 


(2 


where  Ws(f)  is  the  spectral  density  of  the  signal  as  shaped  by  the  receiver. 


Te  =  rms  integration  time 


1 


(2-23) 


T 


2 


1  1 
— + — + 

Yi  Y2 


Y.Y2 


where  —  and  —  are  the 

Y,  Y2 


SNR(s)  for  each  receiver  respectively. 


D.  PERFORMANCE 

Given  (2-20)  and  (2-21),  it  is  clear  that  both  the  GCC  and  the  CAP  can  be 
rendered  ineffective  by  moderate  in-band  noise.  Figure  2-4  below  is  an  illustration  of  the 
theoretical  standard  deviation  of  a  TDOA  measurement  made  at  various  SNR(s)  and 
integration  times.  It  considers  only  the  errors  introduced  by  AWGN  according  to  the 
theory  presented  above  and  includes  no  error  from  any  other  source  such  as 
instrumentation  or  machine  precision  for  example.  Because  theory  dictates  these  values 
as  the  minimum  errors,  they  may  be  considered  the  lower  bound  for  this  scenario. 

The  model  for  this  demonstration  assumes  a  rectangular  signal  spectmm  with 
rectangular,  full  bandwidth  Gaussian  noise  with  no  Doppler  between  observers.  It 
assumes  no  magnitude  or  phase  mismatch  between  receivers  and  the  same  SNR  at  both 
receivers.  It  approaches  the  ideal  situation  given  the  perfect  match  between  signal 
bandwidth,  noise  bandwidth  and  receiver  measurements.  In  reality,  the  signal  bandwidth 
would  be  more  Gaussian  or  raised  cosine  pulse  shaped  and  the  receivers’  measurements 
would  not  be  matched  in  phase  or  magnitude.  These  differences  would  corrupt  the 
TDOA  estimation  further,  requiring  more  SNR  and  integration  time  to  achieve  the  lower 
bounds  determined  by  this  idealized  model. 

Clearly,  below  approximately  10  dB  SNR  and  400  ms  the  idealized  model 
exceeds  a  TDOA  standard  deviation  of  100  m.  Given  three  TDOA  measurements  at  this 
level  and  a  simple  Euclidean  norm  combination  of  these  errors,  a  geolocation  derived  in 
the  absence  of  any  geometric  dilution  (see  Chapter  V)  would  have  a  standard  deviation  on 
the  order  of  173  m.  This  exceeds  the  100  m  goal  of  the  ARL:UT  project  without 
considering  the  additional  errors  introduced  in  an  actual  system  (again  see  Chapter  V). 
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2500 


Figure  2-4  -  Theoretical  standard  deviation  of  TDOA  measurement  given  various 
SNR(s)  and  integration  times 

Poor  SNR(s)  and  SNIR(s)  are  obviously  problematic  for  the  GCC  and  CAF  given 
the  need  for  geolocation  accuracy  on  the  order  of  100  m.  While  increased  integration 
times  can  solve  some  low  SNR  problems,  this  solution  has  an  upper  bound.  The 
integration  time  must  not  exceed  the  coherence  of  the  SOL  That  is,  the  statistics  of  the 
SOI  must  remain  relatively  stationary  during  the  integration  time  for  the  GCC  and  CAF  to 
operate  properly.  A  significant  increase  in  integration  time  to  combat  poor  SNR  may 
exceed  the  coherence  time  of  the  SOI  and  introduce  yet  more  error.  Modeling  the  SOI  as 
cyclostationary  vice  stationary  and  employing  the  appropriate  processing  techniques  can 
in  many  cases  overcome  most  of  these  problems. 

Cyclostationary  techniques  exploit  periodicities  introduced  to  man-made  signals 
in  a  number  of  ways.  These  periodicities  can  unique  to  specific  signals  and  thus  can  be 
used  to  distinguish  one  signal  from  another  in  the  same  bandwidth  as  well  as  significantly 
reduce  the  level  of  post-processing  noise.  Thus,  with  cyclostationary  signal  processing,  it 
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is  possible  to  tolerate  significantly  lower  SNR(s)  and  still  obtain  excellent  TDOA 
measurements.  Chapter  III  introduces  the  theory  of  cyclostationarity  in  communications 
signals  and  develops  a  method  for  TDOA  determination. 
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III.  CYCLOSTATIONARY  TDOA  DETERMINATION 


A.  DEFINITIONS 

The  concepts  of  cyclostationarity  have  been  examined  in  theory  for  over  two 
decades.  Beginning  in  the  late  1960’s,  Dr.  Franks  of  the  University  of  Massachusetts  at 
Amhurst  and  Dr.  Gardner  of  the  University  of  California  at  Davis  (UCD)  began  extensive 
research  in  the  area  of  cyclostationary  signal  processing.  Dr.  Gardner  has  since  produced 
multiple  publications  in  the  field.  His  paper  on  the  general  theory  of  cyclostationary 
signal  processing,  published  in  the  April  1991  edition  of  the  IEEE  Signal  Processing 
Magazine  [6],  stands  as  the  key  reference  for  most  aspects  of  the  theory.  As  is  always  the 
case  when  one  individual  has  contributed  so  much  to  an  important  engineering 
development,  much  of  the  theory  and  examples  presented  here  follow  Dr.  Gardner’s  work 
closely.  His  solo  efforts  and  collaboration  with  others  appear  in  references  [6],  [7],  [8] 
and  [9]  and  is  redundant  in  many  instances.  Where  appropriate,  the  specific  source  of 
unique  information  is  referenced  below. 

As  previously  noted,  most  modem  signal  processing  techniques  associated  with 
communications  applications  treat  SOI(s)  as  stationary  random  processes.  However, 
because  most  manmade  signals  are  generated  through  some  repetitive,  periodic  process 
such  as  the  amplitude,  frequency  or  phase  modulation  of  a  sinusoidal  carrier,  the 
encoding  of  data  or  the  encryption  of  a  message,  their  statistics  inevitably  vary 
periodically  with  time.  While  in  many  instances,  receivers  may  successfully  ignore  the 
underlying  periodicity  of  a  manmade  signal,  often  the  detection  of  the  signal  and  the 
estimation  of  its  parameters  is  more  successfully  accomplished  by  modeling  the  signal  as 
cyclostationary  vice  stationary. 

Simply  stated,  a  process  with  statistics  that  vary  periodically  with  time  is  termed 
cyclostationary.  Figure  3-1  depicts  a  block  diagram  of  the  procedure  that  leads  to  a 
cyclostationary  signal  for  most  communications  processes.  A  stationary  random  message 
such  as  digital  data  or  analog  voice  is  modulated,  clocked  or  framed  by  a  periodic 
communications  process  through  some  nonlinear  system  and  a  cyclostationary  signal 
results. 
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Stationary 
Random  Message 

-  digital  data 

-  thermal  noise 

-  analog  voice 

Cyclostationary 
Signal 

Source  of  Eeriodicity 

-  RF  oscillator 

-  bit-stream  clock 

-  repeating  frame  structure 

-  rotating  machinery  waveforms 

(motors,  engines,  turbines,  propellers...) 


Figure  3-1  -  Origins  of  cyclostationarity,  adapted  from  [10] 

From  a  strictly  mathematical  point  of  view,  a  cyclostationary  signal  of  order  n  is 
one  that  will  have  additive  sine-wave  components  that  result  in  spectral  lines  for  some  n* 
order  nonlinear  transformation  of  the  signal.  In  the  case  of  n  =  2,  a  signal  is  said  to  be 
second  order  cyclostationary  if  a  quadratic  transformation  produces  additive  sine-wave 
components  that  generate  spectral  lines.  This  characteristic  may  be  thought  of  as  akin  to 
a  process  being  considered  wide  sense  stationary  or  stationary  through  order  2. 

To  continue  the  case  of  n  =  2  more  specifically,  a  signal  x(t)  is  cyclostationary 
with  cycle  frequency  a  if  and  only  if  some  of  its  delay  products  y(t)  =  x(t)x(t  -  x)  produce 
a  spectral  line  at  frequency  a.  If  not  all  cycle  frequencies  a  for  which  x{t)  exhibits 
cyclostationarity  are  harmonics  of  a  single  fundamental  frequency,  then  x(t)  is 
polycyclostationary.  Polycyclostationarity  implies  the  existence  of  more  than  one 
periodicity  in  the  statistics  of  a  signal.  This  in  turn  implies  more  than  one  source  of 
periodicity  such  as  the  case  when  clocked-digital-data  phase  modulates  a  sinusoidal 
carrier  in  an  M-ary  phase  shift  keying  scheme  [6]. 
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To  further  illustrate  the  concept,  consider  a  signal  x{t)  that  contains  an  additive 
sine-wave  component  at  frequency  a  and  is  of  the  form 

flcos(2ar  +  6)  with  a  9^0  (3-1) 

The  Fourier  coefficient  defined  as 


M“  = 


(3-2) 


where  (•)  denotes  time-averaging,  will  be  non-zero.  Note  that  this  is  of  the  form  of  the 
common  representation  of  the  Power  Spectral  Density  (PSD)  of  x{t)  with  a  spectral  line  at 
+a  and  its  mirror  at  -a.  In  simple  terms,  x{t)  contains  first-order  periodicity  at  frequency 
a  given  by 


K“r[^(/-«)+^(/+«)]  (3-3) 

Now,  considering  contributions  to  x{t)  from  other  sources,  the  total  signal  may  be 
represented  as 


x(0  =  «cos(27car-t-0)+n(r)  .  (3-4) 

where  n{t)  can  be  thought  of  as  the  random  energy  outside  that  of  the  SOI.  If  n{t)  is 
strong  in  comparison  to  the  SOI  such  that  it  masks  the  sine-wave  components  in  x{t)  from 
detection  during  casual  inspection  of  the  waveform,  then  x{f)  can  be  thought  of  as 
possessing  hidden  periodicity.  This  hidden  periodicity  can  still  be  exploited  by  using 
signal  processing  techniques  such  as  the  PSD  function  as  noted  above.  However  more 
powerful  cyclostationary  signal  processing  techniques  are  more  sophisticated  and  can 
unmask  periodicities  in  signals  that  are  hidden  even  from  common  traditional  techniques 
like  the  PSD  function. 
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B.  CYCLIC  AUTOCORRELATION  FUNCTION 


Traditionally,  a  common  second  order,  time-domain  statistic  used  in  signal 
processing  is  the  autocorrelation  function  given  by  the  quadratic  transformation 

R^{t)  =  {x(t)x(t-r))  (3-5) 

In  the  case  of  a  signal  that  can  be  modeled  as  cyclostationary,  this  transform  will  produce 
spectral  lines  at  non-zero  frequencies  a  such  that 

(3-6) 

where 

y{t)  =  x{t)x{t-x)  (3-7) 

This  signal  x{t)  as  noted  above  contains  second  order  periodicity  manifested  in  the 
PSD  of  the  delay  product  given  in  equation  (3-5)  above.  Transforming  (3-5)  into  a 
symmetric  delay  product  and  accommodating  complex  signals  as  well  gives 

(3-8) 

which  forms  the  basis  for  the  fundamental  second  order  cyclic  moment  known  as  the 
cyclic  autocorrelation  function: 

(3-9) 

Of  note,  the  cyclic  autocorrelation  function  assumes  the  form  of  the  Fourier 
coefficients  of  the  additive  sine-wave  components  produced  by  the  periodicity  of  the 
delay  product  in  (3-6). 
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Another  interpretation  of  (3-9)  is  as  the  traditional  stationary  autocorrelation 
function  multiplied  by  a  kernel  that  produces  spectral  lines  at  frequencies  where  the 
stationary  autocorrelation  function  contains  additive  sine-wave  components  indicative  of 
its  periodicity.  Consequently,  for  either  interpretation,  at  a  =  0,  the  cyclic  autocorrelation 
function  reduces  to  the  traditional  autocorrelation  function. 

A  final  interpretation  of  the  cyclic  autocorrelation  function  can  be  seen  by 
defining  first 


u(t)  =  x(t)e  ^ 


(3-10) 


and 


v{t)  =  x(0c^^“ 


(3-11) 


so  that  u{t)  and  v(0  are  frequency  shifted  versions  of  x(t).  Now  /?“(t)  can  be  written 


<(t)  =  (m|  t+-  IV 


t  +  - 
V  2) 


(3-12) 


which  is  the  cross-correlation  of  the  two  versions  of  x{t)  shifted  up  and  down  by 
frequency  oc/2.  In  other  words,  the  cyclic  autocorrelation  function  may  be  viewed  as  the 
correlation  in  the  time  domain  between  two  values  of  x{t)  separated  in  frequency  by  a. 
Consider  as  an  example,  a  Binary  Phase  Shift  Keyed  (BPSK)  signal  given  by 

s(t)  =  A^d(t)cos(2nfj  +  ^)  (3-13) 


where  d{f)  is  the  binary  modulating  wave  form  consisting  of  positive  and  negative 
rectangular  pulses  given  by 
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(3-14) 


d{t)  =  '^d^qit-t^-nT^) 


Approximated  as  a  random  binary  wave,  it  contains  no  first  order  periodicity  and  thus  no 
spectral  lines  in  its  PSD.  Consequently,  the  PSD  of  the  BPSK  signal  s(t)  is  a  scaled  sine 
squared  function  and  will  also  contain  no  spectral  lines  as  is  evident  its  expression 


•^BPSK  (/)  ~  ^c'^b 


^sin^ 


(3-15) 


Thus,  the  autocorrelation  function  of  the  BPSK  signal,  the  inverse  Fourier  transform  of 
(3~15),  will  contain  no  additive  sinusoidal  components  to  indicate  any  periodicity. 
Multiplying  the  conventional  autocorrelation  function  by  the  cyclic  kernel 
cyclic  autocorrelation  function  follows  as 


R? 


27; 


r"  (t)  cosilTif  T)e 


+  ^  ^^y^-j2n[a+2f,]r„  a-2/,  ^^^^-j2n[a-2f,]r„  I 

47;  ‘‘  i 


(3-16) 


where 


=  L  +  f  ~  I  (3-17) 

This  clearly  shows  additive  sine-wave  components  at  a  =  ±2/^  ±  7?^  and  a  =  ±7?^  [9], 
Thus,  the  quadratic  transform  that  is  the  cyclic  autocorrelation  function  unmasks  hidden 
periodicity  in  this  simple  BPSK  signal  and  proves  it  polycyclostationary  in  the  process. 
Naturally,  the  conventional  cross-correlation  function  and  the  cyclic  cross-correlation 
functions  are  related  in  a  similar  manner. 
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c. 


SPECTRAL-CORRELATION  FUNCTION 


The  frequency  domain  equivalent  for  the  cyclic  autocorrelation  function  is  the 
spectral  correlation  function  (SCF).  It  follows  directly  from  the  Fourier  transform  of  (3- 
9)  and  has  the  form 


^  1  W+Af/2  * 

(3-18) 

where 


(t+TH  . 

I  x{u)e  ’  ^“du 
Jt-Tl2 


(3-19) 


'  is  the  finite-time  Fourier  transform  of  xit).  More  commonly  used  in  cyclostationary 
signal  processing,  the  SCF  reduces  to  the  conventional  PSD  at  a  =  0  just  as  the  cyclic 
autocorrelation  function  reduced  to  the  conventional  autocorrelation  function  at  a  =  0. 

Continuing  the  BPSK  example  from  above,  the  SCF  can  be  found  directly  from 
(3-17)  or  by  taking  the  Fourier  transform  of  (3-15).  It  has  the  form 


5r(/)  = 


/  +/o  + 


r  r  OL 


y 


for  a  —  nl\  and 


^-j2)rta+2/^]i„ 


a 

*2 


^-j27r[a~2f^]to 

J  2) 

(3-20) 


for  a  =  ±2/^  '^nlT^  with  Q(f)  being  the  Fourier  transform  of  the  keying  envelope  q(t) 

[8]. 
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The  SCF  derives  its  name  from  the  fact  that  it  is  in  fact  a  correlation  of  the  SOI  in 
the  frequency  domain.  The  SCF  at  frequency /o  and  cyclic  frequency  oto  is  merely  the 
correlation  of  two  values  of  the  signal  in  the  frequency  domain  separated  in  frequency  by 
CKo  find  centered  at  frequency /o.  Figure  3-2  shows  this  process  concept  graphically. 


Figure  3-2  -  Spectral  Correlation  Function  illustrated,  from  [10] 

The  SCF  is  plotted  on  what  is  called  the  bi-frequency  plane.  The  plane  is  defined 
filong  one  axis  as  spectral  frequency  f  and  along  the  opposing  axis  as  cyclic  frequency  oc. 
Figure  3-3  illustrates  this  point.  The  magnitude  of  the  SCF  corresponds  to  a  height  above 
the  bi-frequency  plane  and  is  often  plotted  in  that  fashion.  Note  that  the  values  for  a 
range  from  -/s  to/s  while  values  for  spectral  frequency /naturally  range  from  -fJ2  to  /s/2. 
Because  cl  corresponds  to  the  separation  distance  of  the  correlated  values  in  the  frequency 
domain,  its  range  extends  twice  that  off. 

The  SCF  computation  can  be  highly  complex  and  demanding  on  any  processor. 
For  that  reason,  two  algorithms  were  developed  in  [7]  and  designed  specifically  for 
computation  efficiency.  The  Fast  Fourier  Transform  Accumulation  Method  (FAM)  and 
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Figure  3-4  -  SCF  of  some  common  modulated  signals,  adapted  from  [10] 
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the  Strip  Spectral  Correlation  Analyzer  (SSCA)  both  simplify  and  thereby  accelerate  the 
computation  of  the  SCF.  Because  both  the  SCF  and  the  Cross  Spectral  Correlation 
Function  (CSCF)  are  necessary  for  cyclostationary  TDOA  algorithms,  SSCA  plays  an 
important  role  in  this  work.  It  is  presented  in  more  detail  in  Chapter  VI. 

Finally,  in  addition  to  their  role  in  the  TDOA  computations,  the  SCF  and  CSCF 
can  also  be  used  for  purposes  ranging  from  signal  detection  to  characterization  and 
estimation  of  parameters.  Specifically,  the  modulation  type  of  an  unknown  signal  can  be 
found  using  the  SCF  of  that  unknown  signal  and  comparing  it  to  SCF(s)  of  known 
modulation  type.  Figure  3-4  shows  a  small  library  of  SCF(s)  from  several  more  common 
modulation  schemes. 

D.  CYCLOSTATIONARY  TDOA  SIGNAL  MODEL 

The  signal  model  used  to  develop  cyclostationary  TDOA  algorithms  is  very 
similar  to  that  presented  in  Chapter  II.  Recall  the  signals  received  at  two  spatially 
separated  observers  can  be  given  by  (2-1)  and  (2-2)  and  are  below  for  clarity 

x(t)  =  s(t)  +  n,(t)  (3-21) 

and 

y(t)  =  A-sit-D)  +  n2(t)  (3-22) 

The  difference  between  the  two  models  lies  in  the  temporal  and  spectral  relationships 
between  s(t),  n\{i)  and  njit).  In  Chapter  11,  n\{t)  and  «2(t)  were  assumed  to  have  no 
temporally  and  spectrally  components  that  overlapped  the  SOI,  5(t).  In  the 
cyclostationary  TDOA  signal  model,  ni(0  and  mit)  represent  all  SNOI(s)  and  noise 
present  at  the  respective  observers.  They  may  or  may  not  contain  co-channel  interferers 
in  reality,  however  for  purposes  of  this  work,  they  are  assumed  to  contain  interference 
that  in  fact  spectrally  overlaps  the  SOI.  Figure  3-5  portrays  this  situation  graphically. 
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Figure  3-5  -  Cyclostationary  signal  model  scenario 


Given  this  model,  the  cyclic  auto  and  cross-correlation  functions  in  general  can  be 
written  as 


i?,“(T)  =  /?;(T)-h/?“(T) 

(3-23) 

R^(t)  =  1  Ap  R“  (T) 

(3-24) 

(T)  =  AR:  (t  -  (T) 

(3-25) 

just  as  (2-3)  -  (2-5)  represented  the  conventional  auto  and  cross-correlation  functions.  In 
addition,  as  with  (2-6)  -  (2-8),  the  spectral  correlation  functions  are 

s::(f)  =  s:if)  +  s"if)  (3-26) 
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(3-27) 


5“  (/)  =  (/)  (3-28) 

Though  (3-26)  -  (3-28)  contain  the  TDOA  of  the  signal  D,  this  parameter  is  now 
masked  by  the  spectrally  overlapping  components  of  n\{t)  and  «2(0-  Thus,  any  attempt  to 
estimate  D  with  traditional  methods  that  compute  the  above  equations  at  a  =  0,  like  those 
illustrated  in  Chapter  II,  would  result  in  a  cormpted  value.  However,  if  s{t)  contains 
some  cycle  frequency,  Oo,  that  it  does  not  share  with  any  component  of  n\(t)  and  n2{f), 
then  cyclostationary  techniques  can  discriminate  those  contributions  to  D  made  by  s{t) 
from  those  contributions  of  n\{t)  and  ^2(0  that  would  otherwise  corrupt  the  estimate. 
Essentially  this  eliminates  the  SNOI(s).  Thus  a  reliable  estimate  of  D  is  possible  even  in 
highly  corrupt  environments  provided  a  unique  a  exists  for  the  SOI.  Spectral  Coherence 
Alignment  is  a  cyclostationary  TDOA  algorithm  employing  measurements  of  both  the 
SCF  and  CSCF.  It  determines  the  TDOA(s)  for  the  TDOA  processor  developed  here  and 
appears  in  detail  below. 


E.  SPECTRAL  COHERENCE  ALIGNMENT 

Spectral  Coherence  Alignment  (SPECCOA)  was  developed  and  presented  in  [8] 
using  an  ad  hoc  minimum  least  squares  (MLS)  approach.  Using  equations  (3-23)  -  (3- 
25),  and  the  assumption  that  s{f)  contains  a  unique  a,  it  can  be  seen  that 

/?“(«)  =  CR:{u-D)e-^^'^  (3-29) 

The  cross-correlation  term  for  n\{t)  and  n2(0  is  eliminated  by  the  use  of  a  cycle  frequency 

unique  only  to  the  SOI.  An  estimate  of  the  TDOA,  D  that  minimizes  the  sum  of  the 
squares  of  error  magnitudes  between  the  measured  value  of  the  left  side  of  (3-29),  R" 
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and  the  measurement  of  the  right  side  of  (3-29),  R“  with  D  substituted  yields  the  MLS 
optimized  value  for  the  TDOA. 

Mathematically,  from  [8], 


A  A 

D=argmin 

A 

C,T 


(3-30) 


A 

where  Ca(T)  is  the  estimate  of  the  MLS  function.  It  was  further  be  shown  in  [8],  that  the 
solution  of  the  MLS  problem  is  of  the  form 


where 


D 


A 

=argmax 


C„(T) 


(3-31) 


caT)=|j/?;(M)i?;(M-T)‘jM 


(3-32) 


(3-33) 


where  (3-33)  derives  from  (3-32)  through  Parseval’s  relation.  Gardner  and  Chen  go  on  to 
further  prove  that  (3-33)  does  indeed  peak  at  t=  D. 

SPECCOA  derives  its  name  from  the  fact  that  the  peak  in  (3-33)  occurs  by 
maximizing  the  correlation  in/of  the  two  spectral  correlation  functions  through  the 
alignment  of  the  phases  of  the  respective  functions. 

SPECCOA  proves  more  useful  in  tactical  applications  than  other  cyclic  TDOA 
algorithms  (see  [5]  and  [9])  by  virtue  of  its  ease  of  implementation  and  performance  in 
corrupt  environments  [9].  Though  other  cyclostationary  algorithms  consistently  out 
perform  SPECCOA  [5],  they  do  so  at  the  expense  computational  complexity  and  speed. 
Because  this  work  is  intended  to  demonstrate  the  feasibility  of  implementing  a 
cyclostationary  TDOA  algorithm  in  an  operational  tactical  system,  SPECCOA  was  the 
logical  choice  for  its  efficiency  and  performance.  Chapter  IV  discusses  the  algorithm 
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which  utilizes  TDOA(s)  to  determine  the  geolocation  of  an  emitter.  Chapter  V  presents 
the  ARL:UT  prototype  TDOA  geolocation  system.  Chapter  VI  contains  specific 
implementation  issues  for  SPECCOA  in  the  MATLAB®  environment. 
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IV.  TDOA  GEOLOCATION  CLOSED  FORM  SOLUTION 


A.  BACKGROUND 

Given  the  TDOA(s)  from  pairs  of  receivers,  the  problem  now  lies  in  determining 
the  location  of  the  transmitter  from  the  intersection  of  the  isochrons  generated.  Many 
distinct  methods  appear  from  the  geometric  interpretation  of  Schau  and  Robinson  in  [1 1] 
to  the  iterative  solution  proposed  by  Loomis  in  [12]  and  the  closed  form  solution 
illustrated  by  Ho  and  Chen  in  [13].  While  somewhat  similar,  they  all  take  a  specific 
course  in  determining  the  location  of  the  transmitter  given  the  determined  case.  What 
becomes  more  difficult  for  each  of  these  solutions  is  the  determination  of  the  contribution 
of  errors  in  the  measurement  of  the  TDOA(s)  to  the  final  result,  the  geolocation. 

In  [14],  Rusu  develops  a  closed  form  solution  avoiding  the  initial  guesswork 
involved  in  an  iterative  technique.  Unlike  [13],  his  solution  is  completely  general.  In 
addition,  he  shows  that  the  errors  in  the  measurements  may  be  propagated  through  the 
mathematical  model  in  a  linear  fashion,  simplifying  the  determination  of  uncertainty  in 
the  geolocation.  The  following  development  draws  exclusively  from  his  derivation. 

The  problem  of  determining  the  location  of  an  emitter  in  three  dimensions  given 
the  times  of  arrival  (TOA)  at  four  spatially  separated  receivers  has  often  been  treated  as  a 
TDOA  problem.  That  is,  time  is  treated  as  absolute  time.  However,  as  Rusu  points  outs, 
the  mathematical  model  is  invariant  to  time  translation  and  thus  may  be  treated  as  a  TOA 
problem  by  setting  the  arrival  time  at  any  one  of  the  receivers  as  the  origin  of  the 
problem,  t  =  0.  This  simplification  without  loss  of  generality  facilitates  the  development 
of  linear  error  propagation  from  the  initial  measurements  through  the  mathematical  model 
to  the  final  result  as  noted  above. 

In  the  presence  of  moving  emitters  or  receivers,  both  TDOA  and  FDOA  must  be 
determined  by  one  process  as  shown  in  the  development  of  the  CAP  in  Chapter  11.  The 
case  of  FDOA  can  be  similarly  simplified  by  assigning  a  Doppler  shift  of  zero  to  one  of 
the  receivers  reducing  the  problem  to  one  of  frequencies  of  arrival  (FOA).  Fortunately, 
given  the  independence  of  the  TOA  portion  of  the  model  from  the  frequency  of 
transmission  and  Doppler  shifts,  the  TOA  and  FOA  equations  may  be  treated  separately. 
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Though  TDOA  is  the  emphasis  of  this  thesis,  for  completeness  both  the  TOA  and  the 
FOA  solutions  will  be  developed  here. 

B.  MATHEMATICAL  MODEL 

The  explicit  functional  model  relating  the  observable  X  to  the  dependent  variable 
of  interest  Y,  in  the  TOA  case  is  the  solution  of  the  set  of  scalar  equations  of  the  form 
F(X,Y)  =  0.  F  may  be  referred  to  as  the  implicit  functional  model,  X  the  independent 
variable  and  Y  the  dependent  variable.  In  the  development  of  the  TOA-FOA  problem,  X 
is  a  32-dimensional  real  vector  state  space  of  observations  while  Y  is  an  8-dimensional 
status  vector  state  space.  Figure  4-1  illustrates  the  scenario  in  general. 


Figure  4-1  -  TDOA  model 

The  TOA-FOA  case  assumes  that  the  i-th  receiver  located  at  rj  =  (Xj,  yi,  zO  and 
moving  with  velocity  Vi  =  (Uj,  Vj,  Wj)  detects  a  signal  at  time  t,  with  frequency  m,.  Four 
observers  for  the  determined  3-dimensional  case,  constitute  the  observation  state  space 
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vector  of  32-dimensions.  The  parameters  of  the  transmitter,  its  location  and  velocity 
given  by  r  =  (x,  y,  z)  and  v  =  («,  v,  w)  respectively,  its  transmit  time  t  and  frequency  a 
compose  the  8-dimesional  status  vector  state  space. 

The  fundamental  equation  relating  the  signal  travel  time  to  the  separation  distance 
between  observer  i  and  an  emitter  assuming  the  signal  travels  at  the  speed  of  light  is 


^](Xi  -  xf  +  (y,.  -  yf  +  (Zi  -  z)^  =  c(ti  - 1)  (4- 1 ) 

where  ■^(xi-xf  -  yf  +  {Zi  -zf  is  the  Euclidean  distance  between  Fn  and  r  denoted 
hereafter  by  ||r„  -  r||  and  c  is  the  speed  of  light.  Forming  appropriate  vectors  using  a 
single  equation  to  include  all  four  observers  produces  the  implicit  functional  TOA  model 

F;(X,Y)  =  ||ri-r||-c(r,-0,  i  =  l,2,3,4.  (4-2) 

Simple  algebraic  transformation  leads  to 

[(x^-x)^+iy,-yy+iz,-zf-it.-tf}  =0,  i  =  1,2,3, 4.  (4-3) 


In  vector  form  the  TOA  equations  are  defined  as  Ft  =  (Fi,  F2,  F3,  F4). 

The  Doppler  equations  can  be  similarly  treated.  The  fundamental  relationship 
between  the  Doppler  shift  observed  by  receiver  n  and  the  relative  velocities  of  the  emitter 
and  observer  is 


(4-4) 


where  tOd  is  the  Doppler  shift  observed,  co  is  the  transmitted  frequency,  (Vn  -  v)  is  the 

r.  -r 


relative  velocity  and 


T„-r 


is  the  unit  vector  from  the  observer  to  the  transmitter. 
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Again  forming  vectors  based  upon  this  relationship  yields  the  implicit  FOA  functional 
model 


F;,4(X,Y)  =  m,-m 


l-(v.  _v). 


/  =  1,2,3,4.  (4-5) 


defined  as  Fp  =  (F5,  Fe,  F7,  Fg). 

Combining  the  TOA  and  FOA  models  produces  the  TOA-FOA  implicit  functional 
model  defined  as  F  =  (Fj,  Fp). 

C.  DETERMINATION  OF  THE  TRANSMITTER  STATE  VECTOR 

The  determination  of  the  8-dimensional  transmitter  state  vector  involves  finding 
for  X,  the  parameter  Y  that  satisfies  the  combined  TOA-FOA  implicit  functional  model, 
F,(X,Y)  =  0,  /  =  1,2.. .,8,  where  the  functions  F,  are  defined  in  (4-3)  and  (4-5).  The 
solution  of  such  a  system  of  equations  occurs  in  two  steps.  First,  the  TOA  equations  can 
be  solved  for  Yt=  (r,  f).  The  resulting  solution  can  then  be  used  in  the  FOA  equations  to 
solve  for  the  unknowns  Yf=  (v,  (o). 

The  TOA  equations  are  irrational  in  their  form  in  (4-3)  and  thus  must  be  squared 
in  order  to  be  solved.  This  produces  a  set  of  quadratics  for  which  two  sets  of  solutions 
exist  corresponding  to  the  two  roots  of  the  equations.  In  most  cases,  Rusu  points  out,  the 
two  solutions  to  the  rational  quadratics  are  also  solutions  to  the  original  irrational  set  of 
equations.  Rarely,  most  often  in  the  2-dimensional  case,  only  one  of  the  two  solutions 
also  solves  the  irrational  set  of  equations  and  leads  to  a  unique  solution  to  the  geolocation 
problem. 

Since  each  solution  to  the  TOA  case  also  produces  a  unique  solution  to  the  FOA 
equations,  more  often  than  not,  there  are  two  distinct  solutions  to  the  problem.  This  is  the 
classic  case  of  ambiguity  encountered  with  every  TDOA  solution  proposed  thus  far. 

Rusu  handles  this  problem  by  noting  that  information  outside  the  algorithm  must  resolve 
the  ambiguity.  Prior  knowledge  of  the  general  area  of  the  target’s  position  might 
eliminate  one  of  the  solutions.  Multiple  measurements  may  also  reveal  one  solution 
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converging  in  one  location  while  the  other  solution  diverges  and  produces  seemingly 
unrelated  geolocations. 

To  solve  the  TO  A  equations,  (4-3)  must  be  squared  to  produce 

{x^-xf  +  {Zi-zf  -{t.-tf  =0,  i  =  1,2, 3,4.  (4-6) 

The  solution  to  (4-6)  is  found  by  subtracting  the  equations  two-by-two  to  eliminate  x^,  y^, 
^  and  This  produces  a  set  of  linear  equations  in  r  and  t  that  may  be  written  as 


Ar  =  qr  +  s 


(4-7) 


As  the  equations  can  be  subtracted  in  any  order  to  produce  such  an  elimination, 
there  exists  many  possible  resulting  systems  of  equations.  One  more  obvious  possibility 
used  by  Rusu  is 


Xj  —  X2 

<N 

1 

ZJ-Z2'' 

X2-X3 

3^2 -Jj 

Z2-Z3 

y3-y4 

^3  ^4  J 

( t  '—t  \ 

q= 

h  “  h 

V3  U) 


(4-8) 


(4-9) 


(4-10) 


This  system  makes  it  possible  to  solve  for  r  in  terms  of  t  producing  a  t- 
parametrized  solution  to  the  TOA  equations 

r  =  gr-i-h,  (4-11) 


37 


where 


g  =  A'^q  and  h=A''s.  (4-12) 

The  introduction  of  this  result  into  the  k'^  range  equation  from  (4-6)  yields  a 
quadratic  in  time  t.  This  equation,  with  some  algebraic  manipulation  can  be  written  as 

at'^  +  2bt  +  c  =  0  (4-13) 


where, 


aHWr-l 

(4-14) 

=  g  li  +  g  r,+/, 

(4-15) 

c  =  ||h-rjr-?*  (4-16) 

The  two  roots  that  are  the  solution  to  (4- 1 3)  produce  two  values  for  r  when 
inserted  into  (4-1 1).  Rusu  notes  that  the  roots  are  usually  real  but  in  the  rare  cases  when 
they  are  complex,  experience  has  shown  him  that  the  real  part  of  the  complex  roots 
suffices  as  a  solution. 

Once  the  TOA  equations  provide  the  emitter  position  and  time  of  transmission, 
these  values  produce  the  transmitter’s  velocity  and  transmitting  frequency  from  the  FOA 
equations.  If  the  FOA  equations  are  rewritten  as 

£lo  .v-ft).  — =  -i+iio.v.  i  =  l, 2,3,4,  (4-17) 

^iO  ®  ^iO 
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the  FOA  case  becomes  a  non-homogeneous  linear  system  that  may  be  solved  using 
standards  procedures  like  Gauss-Jordan  elimination  [14]. 

D.  ERROR  PROPAGATION 

Perhaps  the  most  valuable  portion  of  Rusu’s  solution  comes  with  the  error 
propagation.  He  developed  a  linear  method  through  which  errors  in  the  measurement  can 
be  propagated  through  the  model  and  result  in  predictable  errors  in  the  fix.  The  following 
derivation  is  based  solely  on  his  development  of  this  method.  Though  he  develops  both 
the  TOA  and  the  FOA  propagations,  only  the  TOA  case  will  be  presented  here. 

The  key  to  the  linear  propagation  of  errors  is  based  upon  the  assumptions.  First, 

Y  =  Y(X)  is  assumed  to  be  differentiable  in  the  neighborhood  of  the  measured  value  X 
containing  the  exact  (error  free)  value  of  the  observable  X.  Second,  if  the  errors  in  the 
observation  and  the  parameter  are  be  denoted  5X  =  X  -  X  and  5Y  =  Y  -  Y,  Y  is  assumed 
to  be  linear  in  the  region  of  these  errors  and  thus  6Y  can  be  written 


5Y  = 


,ax, 


5X 


(4-18) 


The  following  Jacobian  matrix  represents  the  derivative  of  a  vector  valued 
function,  Y  =  (Yi,Y2,-,Yn)  with  respect  to  a  vector  variable  X  =  (jci,X2,...,Xn) 


UxJ 


a 


y 


(4-19) 


Given  the  linear  approximation  in  (4-18),  the  covariance  matrix  of  the  parameter 
can  be  related  to  the  covariance  matrix  of  the  observation  through 
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(4-20) 


dxj  \dxj 


The  implicit  function  theorem  provides  a  convenient  method  for  computing  the 
derivatives  of  Y  with  respect  to  X  from  the  relation  F(X,Y)  =  0.  This  provides  the 
general  formula 


UxJ"  iaYj  Uxv 


(4-21) 


Just  as  the  relationship  of  the  TOA  and  FOA  solutions  allowed  the  separation  of 
the  TOA  case  from  the  FOA  case,*the  matrices  involved  in  (4-21)  can  be  written  such  that 
■  the  TOA  and  FOA  contributions  are  separated  as 
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(4-23) 


(4-24) 


The  independence  of  the  TOA  solution  from  the  FOA  equations  is  evident  in  the 


zeros  found  above  corresponding  to 
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The  problem  of  inverting  the  Jacobian,  3F/3Y ,  can  be  solved  numerically 
however,  Rusu  managed  to  express  the  inverse  of  dF/dY  in  terms  of  the  inverses  of  the 
TOA-FOA  blocks  forming  it.  He  points  out  that  this  formally  separates  the  TOA  and 
FOA  error  propagations  and  allows  computation  of  just  the  TOA  solution  and  errors  if 
just  the  transmitter  location  is  desired.  The  block  inverse  of  (4-24)  is 
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(4-25) 


It  is  now  possible  to  write  the  components  of  the  right  hand  side  of  (4-22)  as 
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(4-26) 
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using  (4-21),  (4-23)  and  (4-25). 

Now  the  TOA  Jacobian  matrices  related  to  the  TOA  error  propagation  can  be 
explicitly  expressed  using  the  functional  model  given  in  (4-3).  The  Jacobians  can  be 
expressed  as 
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Implementation  of  Dr.  Rusu’s  algorithm  into  MATLAB®  followed  naturally 
given  its  matrix  and  vector  nature.  It  plays  an  integral  part  in  the  ARL:UT  prototype 
TDOA  geolocations  system  which  is  discussed  in  Chapter  V.  Specifics  regarding  the 
encoding  of  Dr.  Rusu’s  algorithm  can  be  found  in  Chapter  VI. 


42 


V.  THE  CARRY-ON  MULTI-PLATFORM  GPS-ASSISTED  TIME 
DIFFERENCE  OF  ARRIVAL  SYSTEM 

A.  PROJECT  BASIS 

Given  the  problem  of  developing  a  portable,  GPS  assisted  TDOA  system  with 
commercially  available  products,  the  project  team  at  ARL:UT  set  out  to  define  the 
problem  and  develop  solutions  based  on  which  portions  of  the  problem  could  be  most 
influenced.  The  descriptions  and  performance  reports  below  are  based  upon  [15]  and 
numerous  conversations  with  team  members. 

The  TDOA  problem  can  be  reduced  to  three  sub-problems:  geometry,  signal 
timing  and  measurement.  All  three  are  related  to  the  rms  error  of  the  geolocation  by 

^position  ~  GDOP 

The  Geometric  Dilution  of  Precision  (GDOP)  is  a  scalar  multiplicative  factor 
which  serves  to  magnify  any  timing  and  measurement  error.  It  is  directly  related  to  the 
relative  positions  of  the  observers  and  the  emitter.  Poor  GDOP(s),  generally  numbers 
greater  than  three  in  this  application,  result  from  geometries  where  the  observers  are  lined 
up  or  clustered  in  one  area  with  respect  to  the  emitter.  Conversely,  good  GDOP(s), 
numbers  smaller  than  about  three,  result  from  situations  where  the  observers  are 
positioned  such  that  their  TDOA(s)  intersect  at  nearly  right  angles.  In  three  dimensions, 
the  most  favorable  GDOP  occurs  by  maximizing  the  volume  of  the  polyhedron  formed  by 
pointing  four  vectors  from  the  four  observers  to  the  emitter  and  enclosing  the  figure  that 
results. 

The  geometry  of  the  situation  is  essentially  out  of  the  control  of  the  TDOA  system 
as  the  location  of  the  emitter  is  not  known  a  priori  but  would  be  needed  in  order  to 
position  the  observers  in  an  optimal  configuration.  Only  when  many  more  observers 
participate  than  the  minimum  needed  for  a  determined  solution  (three  are  needed  for  2D 
and  four  for  3D)  can  the  geometry  of  the  situation  be  better  controlled  by  the  TDOA 
system.  In  practice,  this  occurs  rarely  and  thus,  the  best  course  of  action  considering  (5-1) 
above  is  to  minimize  the  timing  and  measurement  error  to  the  greatest  extent  possible. 
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GPS  provides  the  capability  to  reduce  both  these  error  sources  significantly.  It  was  with 
this  in  mind  that  the  ARL:UT  project  team,  led  by  Mark  Leach,  pursued  the  design  and 
implementation  of  a  new  GPS-assisited  TDOA  system. 

NSGSA  in  its  tasking  to  ARLrUT  added  some  additional  requirements  in  the  area 
of  component  sources,  architectures  and  interfaces  with  existing  equipment.  The 
following  sections  outline  the  system  in  progressively  more  detail  starting  with  a  general 
summary  and  concept  of  operations,  continuing  into  a  hardware  description  and  ending 
with  a  more  detailed  discussion  of  error  sources  and  system  performance  in  operational 
tests. 


B.  PROJECT  OVERVIEW 

The  system  would  receive  only  a  frequency  of  interest  in  the  HF,  VHF  or  UHF 
line  of  sight  communications  frequency  bands.  It  must  utilize  an  open  system 
architecture  and  maximize  the  use  of  commercial  and  government  of  the  shelf  hardware 
(COTS  and  GOTS  respectively).  It  was  permitted  just  a  single  voice  grade  channel  for 
communication  between  observers  and  needed  to  interface  with  existing  Navy  hardware 
for  that  purpose  (i.e.  ANAVSC-3  transceiver).  Finally,  the  system  needed  to  be 
compatible  with  the  Navy’s  Unified  Build  (UB)  environment  and  interact  seamlessly  with 
the  Joint  Maritime  Command  and  Information  System  (JMCIS). 

The  system  in  its  current  configuration  consists  of  a  network  of  observers 
(minimum  of  three  for  a  geolocation  in  2D  and  four  for  3D)  with  one  observer  acting  as 
the  master  node,  the  others  as  slave  nodes.  It  employs  the  distributed  processing 
technique  depicted  in  Figure  5-1.  Once  given  the  frequency  of  interest  by  an  existing 
signal  acquisition  system  ,  the  master  node  tunes  its  receiver  to  that  frequency  and 
notifies  the  slaves  to  do  likewise.  Each  node  samples  the  incoming  signal  and  buffers  the 
data  in  mass  memory  storage.  The  master  logically  determines  the  400  ms  portion  of  its 
samples  that  contains  the  most  energy  from  the  SOI  and  computes  the  Fast  Fourier 
Transform  (FFT)  of  this  400ms.  The  Master  then  broadcasts  the  FFT  to  the  slaves  with  a 
GPS  time  stamp. 
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Remote  1 


Figure  5-1  -  Distributed  processing,  from  [15] 


Once  the  slaves  receive  the  FFT  from  the  master  they  search  their  buffers  for  the 
400  ms  of  samples  that  correspond  to  the  epoch  of  the  time  stamped  FFT.  They  then  take 
the  FFT  of  their  data  and  perform  a  C AF  to  extract  the  TDOA.  In  addition  to  this  TDOA, 
the  slaves  send  their  GPS  derived  position,  time  and  velocity  (PVT)  measurements 
corresponding  to  the  time  of  intercept  to  the  master. 

The  master  collects  the  TDOA(s)  and  PVT(s)  and  forwards  the  package  to  a 
standard  Navy  TAC3/4  workstation  which  calculates  the  geolocation. 


C.  HARDWARE  DESCRIPTION 

The  system  is  comprised  of  two  major  sub-systems:  the  measurement  subsystem 
and  the  geolocation  subsystem.  Each  is  designed  in  modular  fashion  accepting  certain 
formatted  inputs  and  producing  certain  formatted  outputs.  This  modularity  facilitates 
changes  internal  to  either  system  including  entire  components  or  software  modules  as  the 
process  within  each  of  these  subsystems  is  inconsequential  to  the  other  module.  Only  the 
inputs  and  outputs  are  relevant. 

The  measurement  subsystem  includes  the  hardware  and  software  used  to  compute 
the  TDOA(s).  It  is  based  on  a  13  slot,  C-size  VXI  chassis  housing  an  HP  Signal  Analyzer 
package.  This  package  includes  a  Motorola  68040  based  VXI  controller,  a  10  MHz,  23- 
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bit  (18  linear)  analog  to  digital  converter  (A/D),  a  digital  signal  processing  (DSP)  card, 
200  MB  of  disk  storage  and  software.  The  DSP  module  is  a  Motorola  96002  based 
floating  point  DSP  chipset  on  a  VXI  card.  Also  on  the  chassis  are  a  Ball-Efratom  model 
FRK  10  MHz  Rubidium  (Rb)  reference  oscillator,  a  Precise  Positioning  Service  (PPS) 
capable  GOTS  GPS  receiver  by  Ashtech  Inc.  and  E-systems  and  an  ARL:UT  developed 
clock  synchronization  module.  The  phase  lock  loop  (PPL)  module  allows  the  Rb 
reference  oscillator  to  be  in  synch  with  the  GPS  1  pulse  per  second  (pps)  timing  output. 
The  bulk  of  this  subsystem  is  shown  in  Figure  5-2. 


Figure  5-2  -  Hardware  configuration,  from  [15] 

The  geolocation  subsystem  includes  a  TAC3/4  computer  on  which  reside  the 
algorithms  both  to  compute  the  geolocation  and  communicate  that  fix  to  the  JMCIS.  In 
addition,  it  also  holds  the  system’s  capability  for  simulation  and  playback  modes.  The 
TDOA  algorithms  consist  of  the  closed  form  solution  presented  in  Chapter  V  and  a 
Kalman  filter  method.  The  closed  form  or  determined  solution  is  the  primary  solution 
and  lends  itself  quite  naturally  to  situations  where  the  emitter  has  short  transmission 
durations  which  preclude  the  use  of  the  Kalman  filter  method.  In  the  event  of  long 
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duration  transmission,  the  determined  solution  will  feed  the  Kalman  filter  a  first 
prediction  and  allow  the  Kalman  filter  to  produce  fixes  as  a  function  of  time  over  the 
duration  of  the  transmission.  The  Kalman  filter,  when  given  sufficient  data  and  a  starting 
point,  has  better  error  prediction  capabilities  and  is  thus  the  preferred  method  in  the 
correct  situations. 

The  interface  with  JMCIS  provides  the  system  a  data  display  capability.  By 
updating  the  track  database  through  a  standard  UB  message,  the  system  is  able  to  both 
communicate  with  the  Navy’s  latest  command  and  information  system  and  avoid  having 
to  produce  its  own  man-machine  interface  on  this  output  level. 

Finally,  this  system  is  also  capable  of  simulating  scenarios  in  order  to  investigate 
the  effects  of  variations  on  the  many  factors  involved  from  the  receipt  of  the  signals  to  the 
processing  of  the  data  to  the  computation  and  display  of  the  fix.  In  addition,  a  capability 
to  playback  sessions  recorded  and  stored  on  tape  is  built  in  to  the  geolocation  subsystem. 
This  feature  allows  for  post  processing  and  analysis  of  operational  tests  enabling  finer 
detailed  analysis  of  many  of  the  variables  including  error  source  investigation. 

Perhaps  the  one  key  point  that  stands  above  the  other  accomplishments  of  the 
system  is  the  significant  reduction  of  signal  timing  errors  between  the  collection  nodes. 
With  the  use  of  the  PPS-capable  GPS  receivers,  the  1  pps  output  from  the  GPS  satellites 
in  concert  with  differential  GPS  techniques  enables  the  time  offset  between  two  GPS 
receivers,  separated  by  up  to  hundreds  of  kilometers,  to  be  on  the  order  of  20  ns.  In 
addition,  the  PPL  circuitry  allows  each  node’s  Rb  standard  to  be  synchronized  to  within  5 
ns  of  this  1  pps  output  giving  a  total  signal  timing  error  of  a  mere  25  ns.  In  terms  of  light 
distance,  that  translates  to  7.5  meters.  Figure  5-3  shows  the  entire  system’s  functional 
relationships  in  summary. 
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Figure  5-3  -  Tactical  TDOA  Functional  Diagram,  courtesy  ARL:UT 


D.  ERROR  ANALYSIS 

Many  sources  of  errors  exist  in  the  operational  system.  Another  accomplishment 
of  the  project  is  the  detailed  identification  and  prediction  of  the  magnitude  of  those  error 
sources  which  have  been  validated  by  the  operational  tests  outlined  in  the  next  section.  In 
their  investigations,  the  ARL:UT  team  was  able  to  collect  the  sources  into  four  major 
categories:  fundamental  limits,  instrumentation  limits,  operational  environment  and 
network  geometry. 

The  lower  bound  on  timing  and  measurement  accuracy  may  be  related  to  the 
variance  of  the  TDOA  estimator.  This  in  turn  is  related  to  the  incoming  SNR,  signal  and 
noise  bandwidth  and  integration  time  as  first  shown  by  Hertz  and  Azaria  in  [2]  and  later 
extended  more  specifically  to  this  system  by  Rusu  and  Giulianelli  in  [16]. 

Instrumentation  limits  add  significantly  to  the  errors  in  signal  arrival  time 
determination.  Specifically,  though  ideally  each  GPS  receiver  will  output  the  GPS  1  pps 
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at  the  same  time  during  each  one  second  epoch,  there  will  be  some  offset  in  reality.  In 
addition,  that  1  pps  timing  pulse  is  used  to  trigger  the  sampling  of  the  incoming  signal. 
Delays  from  the  time  the  A/D  converter  receives  the  1  pps  and  the  onset  of  sampling  for 
that  clock  cycle  will  vary  from  one  A/D  to  another.  This  compounds  the  1  pps  offset 
error  and  translates  directly  to  error  in  the  cross-correlation  computation  of  the  TDOA(s). 
The  cross-correlation  methods  of  TDOA  determination  measure  the  time  difference  by 
determining  a  peak  in  the  computed  function.  Five  unknowns  exists  in  this  computation: 
signal  frequency,  signal  bandwidth,  noise  bandwidth,  SNR,  and  statistical  properties  of 
both  signal  and  noise.  The  width  of  this  peak  and  the  presence  of  multiple  peaks  in 
indicative  of  the  magnitude  of  the  lower  bound  on  the  error  in  the  measurement  discussed 
above  and  is  the  direct  result  of  a  combination  of  these  factors. 

The  operational  environment  and  network  geometry  constituted  the  bulk  of  the 
uncontrollable  error  contributions.  While  during  testing,  the  environment  is  somewhat 
controlled,  the  operational  system  would  certainly  be  subject  to  the  conditions  of  the 
current  tactical  situation.  This  holds  true  as  well  for  the  geometries. 

E.  PRELIMINARY  TEST  RESULTS 

Published  test  results  include  both  static  and  dynamic  emitters  as  targets.  In 
November  1994,  the  geolocation  of  a  static  emitter  in  2  dimensions  was  demonstrated  in 
the  Austin,  Texas  area  using  three  observers  with  baselines  of  approximately  33  km.  The 
target  was  a  stationary  NOAA  FM  transmitter  at  162.4  MHz.  The  system  at  this  point 
tested  with  a  2-D  rms  error  of  124  m.  Details  can  be  seen  in  [15]. 

Dynamic  tests  were  conducted  in  January  1995  using  a  moving  emitter,  two 
stationary  observers  and  one  moving  observer.  Again  the  tests  were  conducted  in  the 
Austin  area  and  in  just  two  dimensions  only  given  just  three  observers.  The  target  emitter 
was  an  aircraft  at  approximately  2000  feet,  traveling  at  speeds  of  approximately  100  knots 
and  transmitting  with  a  20  W  UHF  signal  at  461.1  MHz.  The  moving  observer  was  a  van 
mounted  node  traveling  at  speeds  of  up  to  35  mph  in  a  2  square  mile  area.  During  this 
test,  the  prototype  system  performed  geolocations  to  within  277  m  rms  error  using  the 
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Kalman  filter  solution.  This  accuracy  was  achieved  despite  various  geometries  that 
ranged  from  outstanding  GDOP(s)  of  almost  1  to  very  poor  GDOP(s)  of  10. 

While  these  tests  failed  to  achieve  the  goal  of  100  m  rms  error,  instrumentation 
changes  promised  to  increase  the  accuracy  without  significant  changes  to  the  system. 

Two  additional  tests,  each  more  extensive  and  demanding  than  the  previous,  have  been 
conducted  with  the  system.  While  no  performance  results  have  been  formally  published, 
shipboard  tests  off  the  coast  of  San  Diego  California  have  proven  the  system  capable  of 
1 20  m  rms  error  against  representative  VHF/UHF  dynamic  emitters  with  excellent 
geometry  and  -400  m  rms  error  with  poorer  geometries  [17].  Further  changes  in 
receivers  promise  to  remove  biases  in  the  measurements  created  by  the  current  radios  and 
improve  performance  proportionally. 

In  addition,  the  use  of  a  cyclostationary  TDOA  determination  algorithm  should 
provide  some  improvement  in  the  performance  of  the  system  low  SNR  environments. 

The  initial  step  for  the  incorporation  of  cyclostationary  signal  processing  into  the 
ARL:UT  system  is  to  code  the  algorithms  in  MATLAB®  and  test  them  with  ARL:UT  test 
data.  The  coding  of  the  algorithms  in  MATLAB®  appears  in  Chapter  VI. 


VI.  ALGORITHMS 


A.  BACKGROUND 

Any  implementation  of  a  TDOA  processor  utilizing  SPECCOA  must  provide  both 
an  SCF  as  well  as  a  CSCF  at  a  cycle  frequency  Oo,  characteristic  of  the  SOL  In  [19], 
Loomis  and  Berstein  develop  a  functional  TDOA  processor  using  SPECCOA  that 
provides  for  the  efficient  computation  of  both  spectral  correlation  functions  as  well  as  the 
selection  Oq.  This  system  represents  the  basis  for  the  implementation  of  SPECCOA  and 
the  geolocation  algorithm  in  MATLAB®  for  this  project  and  appears  in  Figure  6-1. 


Figure  6-1  -  TDOA  geolocation  processor  functional  diagram,  adapted  from  [19] 

The  four  blocks  of  the  TDOA  processor  constitute  progressively  lesser  degrees  of 
computational  complexity.  The  most  complex  of  the  operations  is  the  calculation  of  the 
SCF(s)  in  the  first  operation  of  the  processor.  The  SCF(s)  of  the  master  and  the  four 
slaves  are  necessary  for  the  selection  of  the  appropriate  cycle  frequency  in  the  next 
operation  of  the  processor,  the  characterization.  The  SSCA  algorithm  computes  the  four 
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SCF(s)  and  uses  a  plotting  routine  to  display  them.  In  the  next  operation,  the 
characterization  is  performed  manually  using  a  plot  of  the  maximum  values  of  the  SCF(s) 
for  each  value  of  cycle  frequency.  This  plot  facilitates  the  choice  of  the  optimum  cycle 
frequency  for  use  in  SPECCOA.  The  third  operation  involves  the  computation  of  the 
SCF  and  CSCF  at  the  cycle  frequency  selected  during  characterization.  Computed  with  a 
simple  frequency  smoothing  algorithm,  the  SCF  and  CSCF  at  the  cycle  frequency  of 
interest  feed  the  final  TDOA  processor  operation,  SPECCOA.  SPECCOA  in  turn 
determines  the  TDOA  by  plotting  the  solution  to  the  ad  hoc  MLS  optimization  found  in 
(3-32).  The  maximum  value  of  this  function  represents  the  TDOA. 

The  SCF(s)  in  the  first  operation  as  well  as  the  characterization  must  be  computed 
for  all  four  observers  for  each  set  of  data.  Cycle  frequencies  for  the  same  emitter  vary 
with  observer  due  to  instrumentation  mismatches  at  each  location.  Details  of  cycle 
frequency  selection  in  these  situations  appear  in  Chapters  VII  and  VIII.  The  cycle 
frequency  specific  SCF  and  CSCF  are  computed  just  three  times  because  as  noted  in 
Chapter  V,  the  only  TDOA(s)  used  in  practice  are  those  between  the  master  node  and  the 
three  slaves.  No  cross  slave  TDOA(s)  are  computed  or  used  in  the  prototype  system  and 
thus  they  are  not  computed  here.  Finally,  three  runs  of  SPECCOA  are  required  to  feed 
the  geolocation  algorithm. 

B.  STRIP  SPECTRAL  CORRELATION  ANALYZER  IN  MATLAB® 

The  SCF(s)  computed  in  the  first  operation  of  the  TDOA  processor  can  place  a 
significant  burden  on  the  processor  if  not  optimized  for  computational  efficiency.  In 
order  to  achieve  the  efficiency  necessary  for  use  in  a  tactical  environment,  a  highly 
expeditious  method  for  computing  and  plotting  the  SCF  is  used.  Directly  from  its 
definition,  (3-17),  the  SCF  can  be  computed  using  a  traditional  time-smoothing  approach 
as  depicted  in  Figure  6-2.  While  this  approach  provides  the  most  accurate  estimate  of  the 
SCF,  it  proves  prohibitively  time  consuming  and  inappropriate  for  most  applications  [20]. 
The  SSCA  algorithm,  a  modification  of  the  time-smoothing  method,  provides  the  best 
combination  of  accurate  estimation  of  cycle  frequencies  of  interest  and  computational 
efficiency  [21].  By  eliminating  one  of  the  band-pass  filters  in  Figure  6-2  and  using  a 
Fourier  transformer  at  the  input  and  output,  SSCA  computes  the  SCF  along 
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Figure  6-2  -  Time-smoothing  spectral  correlation  analyzer  [20] 

diagonal  strips  of  the  bi-frequency  plane  as  depicted  in  Figure  6-3.  While  the  SSCA  does 
manage  to  compute  the  SCF  in  the  most  efficient  manner  developed  thus  far,  it  has  some 
compromise.  The  output  signal  to  noise  ratio  suffers  slightly  as  a  result  of  the  efficiency 
gained  by  the  elimination  of  the  filter  and  the  use  of  the  Fourier  transform.  However,  the 
computational  savings  gained  far  outweigh  the  minor  degradation  in  signal  power  to 
noise  power  at  the  output  of  the  analyzer. 

Figure  6-4  illustrates  the  SSCA  architecture.  The  input  filter  consists  of  a  sliding, 
Hamming-windowed  coarse  FFT  of  length  N' .  Each  of  these  N'  length  segments  is 
essentially  the  output  of  a  band-pass  filter  with  a  bandwidth  of  \/ N' .  This  band-pass 
output  is  then  downcon verted  to  baseband  by  a  complex  exponential  multiplication. 
These  products  are  termed  the  complex  demodulates  of  the  input  signal  and  appear  in 
greater 
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Figure  6-3  -  SSCA  computation  in  the  bi-frequency  plane,  N'  =S 


x(mL  +  n) 


Figure  6-4  -  SSCA  architecture,  from  [20] 


detail  in  Figure  6-5.  As  can  be  seen,  the  input  sequence  is  essentially  decimated  by  the 
process  of  windowing  and  sliding  at  the  input  filter.  Loomis  and  Brown  found  that  the 
most  appropriate  value  for  the  decimation  factor  or  subsampling  parameter  L  was  N' !  A 
[20]. 


Figure  6-5  -  Detail  of  the  computation  of  complex  demodulates 


Because  these  complex  demodulates  are  multiplied  by  the  original  sequence 
before  the  final  FFT,  they  must  be  interpolated  to  the  original  sampling  rate.  In  [20]  and 
[21]  Loomis  and  Brown  use  a  simple  hold  operation.  Linear  interpolation  induces  less 
high  frequency  error  and  adds  little  in  computational  complexity.  It  is  used  here  in  place 
of  the  first  order  sample  and  hold  interpolation.  Finally,  the  interpolated  complex 
demodulates  and  the  original  sequence  are  complex  multiplied  and  the  product  Fourier 
transformed  by  a  full  length,  unwindowed  FFT.  The  relationship  of  the  output  to  the  bi¬ 
frequency  plane  appears  in  Figure  6-6.  Each  row  of  output  represents  one  of 
the  N'  diagonal  strips  in  the  bi-frequency  plane  as  depicted  on  the  previous  page  in  Figure 
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6-3.  Each  column  represents  the  values  of  the  SCF  along  the  diagonal  lines.  The 
resolution  along  lines  of/is  coarse  at  1  /  N'  as  a  byproduct  of  the  computational  savings. 
However,  the  resolution  in  a,  at  1/N,  is  much  higher  which  accommodates  the 
characterization  nicely. 


Figure  6-6  -  Block  form  of  SSCA  output 

Given  that  the  majority  of  operations  in  the  SSCA  can  be  implemented  as  vector 
dr  matrix  operations,  its  implementation  in  MATLAB®  is  fairly  simple.  The  psuedocode 
for  the  algorithm  appears  in  Figure  6-7.  The  MATLAB®  code  for  the  computation  of  the 
SCA  over  the  entire  bi-frequency  plane,  called  ssca,  appears  in  Appendix  A.  Of  note, 
several  MATLAB®  functions  forced  the  use  of  additional  statements  that  might  not 
otherwise  be  found  given  the  psuedocode.  The  fft  function  in  MATLAB®  computes  both 
positive  and  negative  frequencies  but  does  not  place  the  origin  in  the  center  of  the  output 
vector.  Instead,  the  output  offfi  takes  advantage  of  the  cyclic  nature  of  the  DFT.  Its 
output  begins  at  the  origin  of  the  frequency,  moves  through  fs/2  and  then  mirrors  the 
negative  frequencies  with  the  positive  continuing  at  -fs/2  and  ending  with  zero  again. 
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This  forces  any  application  needing  a  double  sided  spectrum  to  rearrange  the  vector  so 
that  the  negative  frequencies  start  the  output  rather  than  end  it.  MATLAB®  has  an 
jftshift  function  that  performs  such  an  operation  and  appears  in  the  code  in  several  places 
for  that  reason.  The  output  of  ssca  is  a  matrix  with  N'  rows  corresponding  to 
the  N'  diagonal  strips  of  the  computation  and  columns  corresponding  to  the  SCF  along 
the  respective  strips.  Referring  again  to  Figures  6-3  and  6-5,  these  strips  are  not 
computed  along  lines  of/ and  a  and  must  therefore  be  transposed  along  those  lines  before 
plotting. 


/*  Compute  Complex  Demodulates  of  x  */ 

Dop:=  0  to  P-1 

Compute  x-j.(pL,4)  =  FFT[a(r)x(pL+r)] 
Dok:=-A^72to  A^72-l 
Compute  XjCpL,  4)  =  xT(pL,  4)exp{-j27tpLk/N' } 

end 

end 

/*  Interpolate  Xt(pL,  4)  */ 

Compute  X^CpL,  4)  =>  (n ,  /^ ) 

/*  Compute  and  Smooth  Product  Waveforms 
Dok:=  -N'flto  N'/2-l 
Compute  S^^(n,f^)  =  y\n)Xj(n,f^) 

Compute  •SJ;(«.A)^  =  FFr{«(n)S5,  («,/,)} 

end 


Figure  6-7  -  Psuedocode  for  SSCA  computation,  from  [21] 

The  program  plotscf  performs  that  function  transposing  the  matrix  output  of  the 
ssca  program  and  plotting  both  the  3-D  SCF  and  the  2-D  magnitude  of  SCF  versus  a  for 
the  characterization  operation.  The  transformation  involves  applying  the  relationship 
between  the  diagonal  strip  subscript  k,  the  frequency /and  the  cyclic  frequency  a.  Those 
relationships  from  [22]  are 
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(6-1) 


N'J 


2  2 


a  =  2f,-2f 


(6-2) 


Given  these  relations  and  the  MATLAB®  5m// function,  the  SCF  plots  easily 
along  lines  off  and  a.  The  2-D  plot  for  additional  characterization  aid  follows  directly 
from  the  SCF  plot. 


C.  SPECTRAL  COHERENCE  ALIGNMENT  IN  MATLAB® 

SPECCOA  represents  yet  another  algorithm  highly  suited  for  MATLAB®.  The 
majority  of  functions  associated  with  computing  (3-32)  are  vector  operations  which 
simplifies  the  MATLAB  code  considerably.  Figure  6-8  is  a  basic  block  diagram  of 
SPECCOA. 


Figure  6-8  -  SPECCOA  block  diagram 
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Once  the  characterization  has  determined  the  cycle  frequency  of  interest,  oto,  a 
frequency  smoothing  algorithm  estimates  the  SCF  and  CSCF  for  oto-  The  inverse  Fourier 
transform  of  the  complex  product  is  then  multiplied  by  the  kernel  .  The  maximum 
value  of  the  real  portion  of  this  product  represents  the  TDOA  of  the  SOI  between  the 
observers. 

SPECCOA  is  used  three  times  to  determine  the  TDOA  between  the  master  node 
and  the  three  slave  nodes.  These  TDOA(s)  along  with  observer  positions  are  inputs  to 
Dr.  Rusu’s  geolocation  algorithm  the  output  of  which  is  the  location  of  the  emitter. 

D.  THE  TDOA  CLOSED  FORM  SOLUTION  IN  MATLAB® 

Dr.  Rusu’s  closed  form  solution  lends  itself  readily  to  MATLAB  implementation. 
The  inputs  include  four  observer  position  vectors  in  WGS84  coordinates  and  three 
TDOA(s),  Di2,  Di3  and  D14.  All  the  intermediate  variables  are  either  vectors  of  length  3 
or  matrices  no  larger  than  3x3.  The  final  solution  involves  finding  the  roots  of  a 
quadratic  that  represent  the  time  of  emission  of  the  SOI  and  substituting  that  time  into  a 
range  equation  to  find  the  WGS84  coordinates  of  the  transmitter. 

Of  particular  note  to  the  MATLAB®  code  that  appears  in  Appendix  A  is  the  use 
of  the  backslash  function,  \,  in  lieu  of  the  matrix  inverse  in  the  algorithm.  This  function 
was  specifically  developed  for  use  in  solving  equations  of  the  form 


Ag  =  q  (6-3) 

for  g  where  A  is  a  matrix  and  g  and  q  vectors.  Normally,  in  matrix  algebra  the  solution 
has  the  form 


g  =  A-'q 


(6-4) 


if  A  has  an  inverse.  However  in  MATLAB®,  it  is  more  accurate  to  use  the  matrix 
division  function  or  backslash.  The  accuracy  of  the  backslash  in  this  scenario  approaches 
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machine  precision  while  the  use  of  the  matrix  inverse  can  reduce  that  accuracy  several 
orders  of  magnitude  [23], 

One  problem  persists  with  the  algorithm  but  is  a  function  of  the  closed  form 
solution  rather  than  its  implementation.  The  two  roots  both  represent  viable  solutions  to 
the  problem.  Information  outside  the  algorithm  must  be  used  to  eliminate  the  ambiguity. 
Multiple  solutions  from  moving  observers  would  more  often  than  not  cause  one  solution 
to  converge  while  the  other  diverged.  No  matter  how  its  resolved  outside  the  algorithm, 
the  ambiguity  always  poses  a  problem  insurmountable  by  the  mathematics. 

The  algorithms  need  to  be  tested  on  controlled  data  prior  to  evaluation  with 
ARL:UT  test  data.  Several  signals  were  generated  using  Simulink®.  These  signals  serve 
to  ensure  the  proper  functioning  of  all  the  code  written  for  this  work.  Following  the 
generated  data,  ARL:UT  data  provides  a  vehicle  for  comparing  the  current  CAF  based 
ARL:UT  method  with  SPECCOA.  The  test  plan  is  detailed  in  Chapter  VII. 


VII.  TEST  PLAN 


A.  GENERATED  SIGNALS 

Three  signals  are  used  in  initially  testing  the  TDOA  processor.  Generated  with 
Simulink®,  all  three  are  BPSK  signals  and  vary  from  noiseless  to  very  poor  SNR.  Each 
signal  was  sampled  at  100  kHz.  They  had  carrier  frequencies  of  0.25  fs  and  data  rates  of 
0.05  fs.  All  were  modulated  with  a  random  binary  signal  generator.  They  are  used  in 
various  combinations  to  verify  the  functions  of  the  three  major  programs,  ssca,  plotscf 
and  speccoa. 

The  first  set  of  tests  involves  using  a  noiseless  BPSK  signal  to  first  verify  ssca 
and  plotscf.  Given  that  these  programs  function  properly,  speccoa  is  tested  finally. 
Figure  7-1  shows  the  simulation  used  to  generate  noiseless  BPSK.  A  total  of  1  second  of 
signal  is  actually  collected  and  stored  as  testbpsk.mat. 


Random  Binary 
Wave 

BPSK  Modulation 

testbpsk.mat 

Figure  7-1  -  Generation  of  noiseless  BPSK  as  testbpsk.mat 
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This  signal  approaches  the  ideal  for  the  processor  for  several  reasons.  First,  it  is 
completely  binary  from  its  inception.  No  analog  to  digital  conversion  is  necessary  and 
thus  introduces  no  error  in  the  process.  Second,  it  is  completely  noiseless.  Finally,  it  is 
modulated  by  a  random  binary  wave  which  introduces  no  biases  in  the  modulated 
waveform.  Given  this,  its  SCF  should  approach  that  given  by  theory  taking  into  account 
the  time  smoothing  effects  of  the  SSCA  algorithm.  Figures  7-2  and  7-3  are  the  time 
domain  and  frequency  domain  plots  respectively  for  testbpsk.mat. 


Figure  7-2  -  Time  domain  plot  of  testbpsk.inat 

The  ideal  SCF  of  a  noiseless  BPSK  signal  can  be  seen  in  Chapter  III,  Figure  3-4. 
The  data  feature  of  testbpsk.mat  should  appear  at  0.05  fs  while  the  twice  carrier  feature 
in  conjunction  with  the  data  feature  appear  at  0.45 /s,  0.5 /s  and  0.55 /s.  Figure  7-4 
displays  the  results  of  ssca  and  plotscf  for  N  =  4096  point  sample.  Only  the  positive 
cycle  frequency  portion  of  the  bi-frequency  plane  is  plotted.  The  cycle  frequency  axis 
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draws  from  upper  right  to  lower  left  and  is  plotted  in  terms  of  nfJN,  n  as  axis  label,  UN  as 
cycle  frequency  resolution.  The  spectral  frequency  follows  from  upper  left  to  lower  right 
labeled  in  the  same  manner.  They  are  unlabeled  in  the  figure  to  eliminate  the  additional 
clutter.  As  noted  in  Chapter  VI,  the  SSCA  algorithm  trades  resolution  in  the  frequency 
domain  for  computational  efficiency  without  losing  the  four  key  features. 


350 


frequency  (Hz) 


x10 


Figure  7-3  -  Frequency  domain  plot  for  testbpsk.mat 

Finally  for  testbpsk.mat,  speccoa  is  run  with  observer  1  collecting  testbpsk.mat 
samples  1  to  8192  and  observer  2  collecting  testbpsk.mat  samples  71  to  8263  for  a  70 
sample  delay.  The  results  of  this  test  appear  in  Figure  7-5.  Because  the  signals  at  the  two 
observers  are  completely  correlated,  this  scenario  represents  an  ideal  situation  for  the 
algorithm.  Its  complete  success  is  expected  for  this  case  but  serves  as  an  initial  check  of 
speccoa.  The  remaining  two  signals  introduce  both  noise  and  poorly  correlated  signals  to 
the  scenario. 
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Two  very  noisy  BPSK  signals  constitute  the  remaining  generated  signal  testing  for 
the  system.  They  were  created  by  passing  the  noiseless  BPSK  signal  through  separate 
AWGN  channels  with  Eb/No  of  0.88  or  -0.534  dB.  This  corresponds  to  probability  of  bit 
error  of  10' .  Figure  7-6  shows  the  generation  of  nlbpsk.mat  and  n2bpsk.inat.  Because 
the  AWGN  channel  introduces  significant  noise  the  BPSK  signals,  nlbpsk.mat  and 
nlbpsk.mat  are  highly  corrupted.  Figures  7-7  and  7-8  show  the  time  domain  and 
frequency  domain  plots  respectively  for  the  two  signals. 


delay  (samples) 

Figure  7-5  -  Output  of  speccoa  routine  for  testbpsk.mat  with  70  sample  delay 

Naturally,  the  level  of  noise  in  these  BPSK  signals  will  be  evident  in  their  SCF(s). 
Figures  7-9  and  7-10  display  the  output  of  ssca  and  plotscf.  Note  that  the  reduction  in 
SNR  compared  to  testbpsk  has  added  noise  to  the  plots  but  has  not  masked  the  important 
features.  While  in  theory  as  presented  in  Chapter  III,  the  AWGN  should  have  no  features 
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other  than  those  evident  at  a  =0,  it  does  add  noise  to  the  higher  values  of  a  in  the  SCF 
plots  in  Figure  7-9.  It  does  so  because  the  AWGN  channel  used  in  the  simulation  is  not 
ideal  AWGN  and  thus  is  not  truly  random  and  Gaussian.  Thus,  some  statistical 
impurities  in  the  channel  as  well  as  the  computational  trade-offs  in  the  estimation 
procedure  of  the  SSCA  algorithm  lead  to  contributions  above  a  =  0. 


Figure  7-6  -  Generation  of  corrupted  BPSK  as  nlbpsk.mat  and  n2bpsk.inat 

Finally,  speccoa  uses  nlbpsk.mat  for  observer  1  and  n2bpsk.mat  delayed  by  70 
samples  for  observer  2  to  verify  that  it  indeed  works  in  highly  corrupt  environments  given 
the  proper  signal  features  as  inputs.  The  results  appear  in  Figure  7-11.  The  successful 
conclusion  of  testing  with  generated  signals  leads  next  to  testing  of  the  ARL:UT  data 
from  a  November  1995  test  in  the  Austin,  Texas  area. 
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Figure  7-1 1  -  Output  of  speccoa  routine  for  nlbpsk.mat  and  n2bpsk.mat 

B.  ARLrUT  COLLECTED  DATA 

The  portion  of  the  November  1995  Austin  tests  that  are  used  in  this  work  were 
collected  by  four  stationary  observers  in  the  Austin  area.  The  emitter  was  a  push-to-talk 
narrow-band  FM  (NBFM),  1  W  transmitter  located  near  the  center  of  the  four  observers. 
The  GDOP  of  the  scenario,  given  the  positions  of  observation  and  the  stationary  target, 
was  excellent  at  less  then  two. 

During  the  test,  each  observer  records  -800  ms  of  data  every  5  seconds.  This  data 
was  stored  to  magnetic  tape  for  post  test  analysis  and  further  off-line  testing.  It  is  used  to 
compare  both  TDOA  determination  and  geolocation  performance  of  the  ARLrUT 
algorithms  and  those  developed  here.  In  order  to  best  understand  the  output  of  ssca, 
plotscf  and  speccoa  for  the  ARL:UT  data,  it  is  necessary  to  follow  the  signals’  processing 
from  antenna  to  storage  media.  The  processing  at  each  observer  is  identical  with  respect 
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to  archiving  the  data  for  later  use  and  thus  a  single  explanation  serves  for  all  four 
observers  in  this  case. 
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Figure  7-12  -  Beginning  of  ARL:UT  signal  processing 

Following  reception,  the  10  kHz  bandwidth  signal  is  heterodyned  to  an 
intermediate  frequency  (IF)  of  21.4  MHz.  At  the  IF,  the  10  Msamples/second  A/D 
converter  aliases  the  signal  at  +1.4  MHz  and  -1.4  MHz  with  additional  spectra  added  at 
intervals  of  the  sampling  rate.  The  spectrum  at  1.4  MHz  is  downconverted  to  baseband 
as  the  product  of  an  appropriate  complex  exponential  multiplication.  Its  mirror  at  -1.5 
MHz  is  shifted  to  -2.8  MHz  in  the  process.  This  series  of  steps  is  summarized  in  Figure 
7-12  with  a  spectrum  relatively  representative  of  that  of  the  test  data. 

At  this  point,  the  data  is  decimated  in  time  by  a  factor  or  1024.  This  step  again  aliases  the 
signal  at  intervals  of  the  new  sampling  rate,/s  =  9765.625  Hz.  The  result  is  a  set  of  four 
spectra  in  the  interval  -/s  /2  to  +/s  /2.  Figure  7-13  shows  the  results  of  this  process  in  the 
time  domain  and  frequency  domain.  Note  the  relationship  between  the  highest  negative 
frequency  spectral  line  at  approximately  -440  Hz  and  the  highest  positive  frequency 
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spectral  line  at  approximately  4440  is  exactly /s/2  as  they  are  related,  aliased  versions  of 
the  spectrum.  The  lowest  negative  frequency  at  approximately  -4460  Hz  and  lowest 
positive  frequency  at  approximately  420  Hz  share  that  same  relationship  for  the  same 
reason.  The  frequencies  are  approximated  above  as  they  change  slightly  throughout  the 
test  and  vary  by  observer  as  well. 

The  corresponding  SCF  has  features  along  the  cyclic  frequency  axis  that 
correspond  to  the  various  combinations  of  separation  between  the  four  spectral  lines  in 
the  PSD  of  the  signal.  The  greatest  contribution  to  the  SCF  occurs  at  a  separation  of/  /2 
given  the  two  constant  relationships  between  the  spectral  lines  sited  above.  At  a  =/ 12, 
all  four  lines  contribute  to  the  SCF  thus  the  a  =/  /2  feature  contains  the  most  energy. 
This  is  best  illustrated  the  SCF  and  the  magnitude  vs.  cycle  frequency  plots  in  Figure  7- 
14.  Note  the  greatest  peak  in  the  SCF  occurs  at/  /2. 

Finally,  using  this  greatest  feature,  speccoa  finds  the  TDOA  between  observers. 
While  other  features  may  also  lead  to  reliable  answers,  the/  /2  feature  is  a  constant 
throughout  the  testing.  All  observers  have  this  feature  and  it  is  always  the  strongest 
feature  of  the  SCF.  An  example  of  the  output  from  speccoa  appears  in  Figure  7-15.  It 
depicts  the  TDOA  between  the  master  node  and  slave  #1.  Specific  results  from  the  testing 
including  some  additional  steps  taken  for  more  resolution  in  the  TDOA  computations 
appear  in  Chapter  DC. 


73 


frequency  (Hz) 


Figure  7-13  -  Time  (top)  and  frequency  (bottom)  plots  of  master  signal,  time  409400 
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Figure  7-15  -  Output  of  speccoa  for  master  and  slave  #3,  time  409510 
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VIII.  GEOLOCATION  SOFTWARE  WORKBENCH  BLOCKS 


A.  SIMULINK®  BLOCK  CONSTRUCTION 

In  addition  to  the  testing  of  the  combination  of  SPECCOA  and  Rusu’s  closed 
form  TDOA  geolocation  solution,  the  algorithms  serve  as  the  first  blocks  readied  for 
construction  for  a  developing  geolocation  software  workbench.  The  details  of  their 
construction  are  outlined  below  in  theory. 

Blocks  that  may  be  used  in  simulations  run  in  the  Simulink®  environment 
conform  to  a  general,  highly  flexible  standard  creates  by  The  Mathworks,  Inc. 
programmers.  They  may  be  implemented  in  one  of  three  ways:  graphically,  through  m- 
files  or  through  another  language  such  as  ANSI  C  which  produces  a  mex-file.  No  matter 
the  method,  the  result  is  an  S-function  that  operates  according  to  the  rules  and  protocols 
of  the  Simulink®  environment. [24] 

Graphical  creation  of  an  S-function  involves  connecting  the  existing  blocks  that 
perform  the  needed  function,  collecting  them  as  a  single  entity  and  masking  that 
collection  as  the  new  block.  This  procedure  forces  Simulink®  to  create  an  S-function 
that  describes  the  new  system.  This  is  by  far  the  easiest  method  of  creation  for  a  speccoa 
block  and  a  closed  form  TDOA  geolocation  block.  The  other  two  methods,  m-files  and 
mex-files,  are  appropriate  in  some  instances  such  as  certain  state-space  systems  or  other 
applications  where  the  system  can  easily  be  written  as  a  set  of  differential  equations. 
Though  more  complex  to  create,  S-functions  created  as  m-files  and  mex-files  simulate 
considerably  faster.  However,  with  the  addition  of  the  Simulink  Accelerator®, 
graphically  created  S-functions  enjoy  some  improvement  in  execution  time  as  well.[24] 

B.  SPECCOA  BLOCK 

All  the  functions  needed  to  use  the  speccoa  code  written  for  the  testing  described 
in  both  Chapters  VII  and  DC,  exist  in  Simulink®  or  one  of  the  toolboxes.  The  Signal 
Processing  Toolbox  ®  and  DSP  Blockset®  contain  the  necessary  buffers  and  functions 
needed  by  the  speccoa  code.  A  user  defined  Matlab®  function  block  serves  as  the 
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vehicle  for  the  speccoa  code  to  be  used  in  the  Simulink®  environment.  Connecting  the 
buffers  with  the  user  defined  function  block  constitutes  the  majority  of  the  speccoa  block. 

In  general,  a  geolocation  workbench  block  of  speccoa  requires  two  scalar 
sequences,  the  time  samples  of  the  two  signals  and  one  scalar  parameter,  the  cycle 
frequency  of  interest  as  inputs.  Time  domain  samples  feed  a  buffer  of  length  defined  by 
the  user  in  order  to  achieve  the  desired  resolution.  These  time  domain  samples  come  to 
the  speccoa  block  from  an  appropriately  designed  data  conversion  block.  Regardless  of 
its  input,  the  conversion  block  must  provide  the  speccoa  block  with  single  real  or 
complex  values  for  each  simulation  time  period.  In  addition,  the  characteristic  cycle 
frequency  enters  the  block  following  the  operation  of  a  characterizer  (defined  elsewhere) 
to  complete  the  inputs.  The  output  of  the  speccoa  block  is  a  2  x  1  vector,  a  time-stamp 
and  a  TDOA.  Figure  8-1  shows  a  schematic  of  the  speccoa  block.  Its  role  in  a  possible 
simulation  scenario  appears  in  Figure  8-3. 


Figure  8-1  -  Schematic  of  speccoa  geolocation  workbench  block 
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c. 


GEOLOCATION  BLOCK 


The  closed  form  TDOA  geolocation  block  is  the  more  complex  of  the  two 
algorithm  blocks.  Its  inputs  include  four  3x1  observer  position  vectors  with  x,  y  and  z 
parameters  for  each  observer.  In  addition,  it  must  have  the  output  of  three  speccoa 
blocks.  The  geolocation  block’s  output  is  a  4  x  1  vector  that  includes  the  target’s  x,  y  and 
z  coordinates  and  a  time  corresponding  to  that  geolocation.  Figure  8-2  illustrates  a  basic 
diagram  of  the  block.  Its  role  in  a  possible  simulation  scenario  appears  in  Figure  8-3  as 
well. 


Figure  8-2  -  Closed  form  TDOA  geolocation  solution  block  schematics 


D.  SIMULATION  SCENARIO 

One  possible  scenario  that  might  be  simulated  using  the  two  blocks  described 
above  is  the  November  1995  Austin  test  of  the  ARLrUT  system.  Using  the  signal  data 
files  and  position  files  for  each  of  the  four  observers,  both  the  speccoa  block  and  the 
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closed  TDOA  form  geolocation  block  are  testable.  The  addition  of  a  characterizer  is 
necessary  to  feed  the  speccoa  block.  This  function  could  be  performed  by  a  user  defined 
block  that  calculates  the  maximum  value  of  the  SCF  for  each  value  of  cycle  frequency  as 
plotscf  plotted  (see  Chapter  VII).  The  maximum  value  of  this  function  could  be  the 
characteristic  cycle  frequency  used  in  the  speccoa  block.  In  addition,  the  geolocation 
block  could  feed  a  plotting  routine  to  display  each  geolocation.  A  diagram  of  this 
implementation  appears  in  Figure  8-3. 


Figure  8-3  -  Simulation  scenario 

Clearly  the  Simulink®  environment  is  ideal  for  testing  the  many  aspects  of 
geolocation  from  data  processing  to  parameter  computation  to  geolocation  estimation. 
The  user  defined  MATLAB®  function  block  provides  for  easy  introduction  of  m-file 
routines  into  the  S-function  domain.  It  allows  the  interchange  of  various  blocks  in  any 
scenario.  This  fact  alone  lends  great  flexibility  and  power  to  a  geolocation  software 
workbench  centered  around  this  environment.  Other  blocks  such  as  pre-filters  for  the 
data  or  logic  to  evaluate  the  quality  of  the  TDOA(s)  or  geolocation  can  also  be 
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constructed.  In  addition,  data  from  many  different  sources  is  handled  simply  by 
modifying  the  data  conversion  block  to  accept  a  different  form  of  input.  Its  output 
remains  the  same.  Using  the  Simulink®  environment  eliminates  the  need  to  define 
standards  for  module  construction  and  provides  an  established  simulation  engine  with 
many  options.  It  is  surely  the  most  expeditious  and  intelligible  implementation  of  a 
geolocation  simulation  system. 
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IX.  RESULTS 


A.  BACKGROUND 

Two  comparisons  are  made  between  the  cyclostationary  TDOA  processor 
measurements  and  those  of  the  ARL:UT  system.  Both  the  TDOA(s)  and  the  resulting 
geolocations  serve  as  benchmarks  of  performance.  In  comparing  the  two  systems’ 
results,  theoretical  TDOA(s)  calculated  from  the  target’s  position  and  the  positions  of  the 
four  observers,  using  a  value  of  2.998  x  10^  for  c,  the  speed  of  propagation,  are  the 
baseline  for  comparison  of  the  TDOA  measurements.  The  target’s  actual  position 
throughout  the  test  acts  as  the  baseline  for  comparison  of  the  geolocations. 

The  TDOA(s)  measured  by  both  systems  are  biased  by  the  radios  used  in  the  test. 
This  constant  bias  was  discovered  in  post  test  analysis  and  removed  from  the  ARL:UT 
published  results.  However,  the  biases  remain  intact  here  to  simplify  the  calculations 
and  comparison.  The  geolocations  suffer  from  the  same  biases  and  are  also  left 
uncorrected  for  the  comparison. 

The  comparison  consists  of  25  time  periods  of  data  from  the  November  1995 
Austin  test.  Because  all  four  observers  recorded  data  during  these  times,  geolocations  in 
three  dimensions  are  possible.  The  process  of  arriving  at  three  TDOA(s)  and  a 
geolocation  involves  plotting  the  SCF(s)  for  the  four  observers,  choosing  a  cycle 
frequency,  feeding  that  data  to  speccoa  and  finally  running  the  closed  form  solution  for 
the  geolocation.  The  plots  for  the  four  SCF(s),  the  four  magnitude  of  SCF  vs  cycle 
frequency  and  the  three  speccoa  outputs  appear  in  Appendix  B  for  each  of  the  25  time 
periods.  As  noted  in  Chapter  VII,  the  feature  of  the  SCF  at /i/2  provides  speccoa  the 
needed  cycle  frequency  in  every  time  instance.  A  summary  of  the  results  appears  below. 

B.  TDOA  COMPARISON 

The  comparison  of  TDOA(s)  between  the  master  node  and  slaves  1,  2  and  3 
appear  below  in  Figures  9-1  through  9-6  below.  Figures  9-1, 9-2  and  9-3  merely  plot  the 
TDOA(s)  in  the  order  in  which  they  were  determined. 
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times 

Figure  9-2  -  TDOA(s)  for  master  node  and  slave  #2,  ARL:UT  (top),  speccoa  (bottom) 
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Figures  9-4,  9-5  and  9-6  show  the  deviation  in  seconds,  from  the  calculated, 
theoretical  TDOA,  for  each  of  the  25  time  periods  along  the  horizontal  axis.  The  top 
plots  in  each  figure  again  show  results  from  the  ARL:UT  test;  the  lower  plots  show 
results  of  the  cyclostationary  TDOA  processor. 

C.  GEOLOCATION  COMPARISON 

While  the  algorithm  in  both  instances  is  identical,  a  comparison  of  the  results 
serves  to  further  illustrate  the  performance  of  the  two  different  TDOA  determination 
algorithms.  The  method  of  comparison  is  merely  the  geolocation  performance. 
Unfortunately,  even  given  the  wi4p  variety  of  TDOA(s)  determined  by  the 
cyclostationary  TDOA  processor,  the  version  of  the  closed  form  TDOA  geolocation 
solution  coded  in  MATLAB®  produced  nearly  the  same  result  for  all  25  time  instances. 
While  this  consistency  is  commendable,  the  various  TDOA(s)  should  have  produced 
fixes  separated  by  several  kilometers.  This  lack  of  definite  variation  in  the  geolocations 
reflects  poorly  on  the  quality  of  the  algorithm’s  outputs. 

The  actual  location  of  the  transmitter  is  given  by  the  coordinates 
Xt  =  -741941.474  m,  yt  =  -5462120.413  m  and  Zt  =  3198242.718  m.  The  MATLAB® 
solution  produces  xm  =  -740591  m,  yM  =  -5443692  m  and  zm  =  3187478  m.  The  only 
change  to  these  rounded  integer  values  over  the  course  of  the  25  time  periods  occurs  in 
the  decimal  places.  Clearly  this  is  a  very  poor  fix  in  addition  to  the  lack  of  variation. 

The  ARL:UT  system  successfully  found  the  coordinates  to  within  100  m  on  all  three  axis 
for  all  25  time  periods.  More  often  than  not,  the  errors  in  the  ARL:UT  calculated  values 
were  on  the  order  of  10  to  20  meters. 

Because  the  MATLAB®  algorithm  did  not  perform  properly,  a  valuable 
comparison  is  impossible.  The  TDOA(s)  discussed  above  are  the  only  viable  measure  of 
comparison  available. 
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Figure  9-5  -  Deviation  from  theoretical  for  TDOA  between  master  node  and  slave  #2, 
ARL;UT  (top),  speccoa  (bottom) 
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D.  COMMENTS 


Clearly  the  cyclostationary  TDOA  processor  does  not  perform  as  predicted  by 
theory  and  as  verified  by  simulation.  Due  to  the  decimation  of  the  original  time  domain 
signals,  several  additional  processing  steps  are  necessary  in  order  to  achieve  the  needed 
resolution  in  the  TDOA  computation.  These  steps  include  resampling  the  signal  at  a 
much  higher  rate  in  order  to  measure  a  TDOA  on  the  order  of  10'^.  This  step  causes 
problems  for  the  speccoa  code. 

The  average  values  for  the  two  TDOA  determination  methods  appear  in  Tables 
9-1  to  9-3  below.  Also  in  the  table  are  the  values  for  the  average  deviation  of  the 
TDOA(s)  from  the  theoretical  values. 

The  code  for  the  geolocation  algorithm  obviously  has  flaws  as  well.  Its  inability 
to  distinguish  between  the  widely  varying  sets  of  TDOA(s)  with  widely  varying 
geolocations  indicates  a  fundamental  programming  flaw. 


Algorithm 

TDOAmi 

Theoretical  TDOAmi 

DeviationMi 

Cyclostationary  TDOA 

-1.308e-5 

1.434e-6 

1.165e-5 

ARL:UT  System 

-1.805e-6 

1.434e-6 

3.66e-7 

Table  9-1  -  TDOAmi  summary 


Algorithm 

TDOAm2 

Theoretical  TDOAmi 

DeviationM2 

Cyclostationary  TDOA 

8.44e-7 

3.828e-6 

2.984e-6 

ARL:UT  System 

2.697e-6 

3.828e-6 

1.131e-6 

Table  9-2  -  TDOAm2  summary 
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Algorithm 

TDOAm3 

Theoretical  TDOAms 

DeviationM3 

Cyclostationary  TDOA 

-5.869e-6 

-3.124e-6 

2.745e-6 

ARL:UT  System 

-3.47  le-6 

3.124e-6 

3.47e-7 

Table  9-3  -  TDOAms  summary 
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X. 


CONCLUSIONS 


A.  OBSERVATIONS 

The  cyclostationary  TDOA  processor  in  combination  with  the  closed  form 
geolocation  solution  in  MATLAB®  does  not  outperform  the  CAP  based  system. 

Because  the  results  deviate  significantly  from  the  theory,  several  additional  tests  are 
necessary  to  isolate  the  cause  of  the  anomaly.  First,  speccoa  must  run  with  fully 
decimated  ARL:UT  data.  Because  the  decimation  exceeds  the  resolution  needed  to 
determine  the  TDOA(s),  properly  the  result  should  be  zero  in  every  instance.  It  is  not. 
The  algorithms  returns  random,  obviously  unreliable  answers  for  the  real  world  test  data. 

The  next  step  is  the  verification  of  zero  delay  with  generated  data.  The 
generated  BPSK  signals  with  zero  delay  produce  zero  as  a  result  with  speccoa.  Thus, 
speccoa  indeed  functions  properly  even  at  delays  below  the  resolution  of  the  sampling 
rate.  However,  it  does  not  function  properly  with  the  ARL:UT  data. 

Finally,  since  the  SPECCOA  algorithm  reduces  to  a  GCC  at  a  =  0  in  the  absence 
of  Doppler  shift,  speccoa  can  be  tested  as  a  GCC  to  reproduce  the  ARL:UT  results.  This 
test  produces  random,  unreliable  results  yet  again  with  ARL:UT  data.  It  does,  however, 
produce  the  correct  results  with  the  generated  BPSK  signals  once  again. 

To  summarize,  plotscf  and  the  cycle  frequencies  it  produces  are  verified  by  the 
correct  plots  of  the  generated  signals  and  the  correct  determination  of  the  dominant 
cyclic  features  of  the  BPSK  signals  even  in  the  presence  of  significant  noise.  Thus,  the 
cycle  frequencies  used  in  speccoa  for  both  the  generated  signals  and  the  ARL;UT  data 
are  deemed  reliable.  Speccoa  on  the  other  hand,  produces  perfect  results  for  the 
generated  signals  at  all  significant  cycle  frequencies  and  at  a  =  0.  However,  it  does  not 
produce  consistent  results  for  the  dominant  cyclic  feature  in  the  ARL:UT  data;  nor  does 
it  produce  results  at  other  cycle  frequencies  that  are  consistent  with  the /s/2  feature’s 
results. 

These  additional  details  lead  to  one  most  likely  conclusion.  The  ARL;UT  data 
does  not  reach  the  algorithms  in  the  correct  form  from  storage  to  processing.  Most 
probably,  the  format  is  incorrect.  In  its  reading  into  the  MATLAB®  environment,  the 
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data  is  essentially  corrupted  producing  such  seemingly  contrary  results.  The  generated 
signals  are  generated  in  the  MATLAB®  environment,  producing  no  such  format 
problems.  The  ARL:IJT  data  is  originally  stored  in  UNIX  format,  byte  reversed  for  PC 
usage  and  then  loaded  into  the  MATLAB®  environment  for  processing.  It  is  this  last 
step  that  is  believed  to  cause  the  anomalies  in  the  results. 

As  a  positive  outcome,  the  testing  did  provide  an  idea  for  the  computational 
complexity  that  a  cyclostationary  TDOA  system  would  require.  The  process  from 
loading  the  signals  to  finding  a  geolocation  required  approximately  2  minutes  per 
geolocation  on  a  Pentium  Pro™  based  machine  with  64  MB  EDO  RAM.  This  is  not 
tactically  feasible.  While  goding  the  system  in  C  would  increase  the  speed  significantly, 
the  computations  of  the  system  would  remain  the  limiting  factor  in  its  implementation. 

The  cyclostationary  TDOA  processor  holds  promise  for  near  future 
implementation.  The  problems  encountered  are  tractable.  All  that  remains  is  the 
computational  power  to  run  the  system  in  a  tactical  environment. 

B.  AREAS  OF  FURTHER  STUDY 

Determining  the  proper  method  with  which  to  load  the  ARL:UT  data  into  the 
MATLAB®  environment  remains  the  highest  priority  for  solving  the  major  problems  of 
the  current  algorithms.  If  this  is  the  cause  of  the  significant  inconsistency  encountered, 
its  solution  will  surely  improve  the  results  to  accuracies  predicted  in  theory. 

The  polynomial  fit  used  by  the  resample  function  in  MATLAB®  needs  to  be 
examined  closely  to  determine  the  extent  of  phase  error  it  introduces  to  the  resampled 
signal.  This  would  allow  for  the  development  of  a  more  accurate  method  of 
interpolating  baek  to  a  sampling  rate  that  provides  sufficient  resolution  for  the  TDOA 
measurement. 

The  closed  form  solution  error  must  be  found  and  corrected.  It  is  likely  that  the 
error  involves  the  reintroduction  of  c  =  2.998  x  10^  once  the  contributions  of  the 
TDOA(s)  are  mapped  to  the  three  coordinates.  In  his  solution,  Rusu  sets  c  =  1  for 
simplicity.  While  this  works  for  the  derivation,  it  most  likely  must  be  given  its  proper 
value  before  the  final  solution  can  be  reached. 
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Once  the  two  algorithms  were  corrected,  using  real  world  signals  from  a  more 
hostile  signal  environment  for  comparison  would  highlight  the  benefits  of 
cyclostationary  processing  best. 

Finally,  further  development  of  Geolocation  Software  Workbench  blocks  would 
accelerate  the  introduction  of  a  simplified  testing  platform  for  all  aspects  of  geolocation. 


95 


96 


LIST  OF  REFERENCES 


1.  Torrieri,  Don  J.,  "Statistical  Theory  of  Passive  Location  Systems",  IEEE  Transaction 
on  Aerospace  and  Electronic  Systems,  Vol.,  AES-20,  Number  2,  March  1984,  p.l83. 

2.  Azaria,  M.,  and  Hertz,  D.,  “Time  Delay  Estimation  by  Generalized  Cross  Correlation 
Methods”,  IEEE  Transactions  on  Acoustics,  Speech  and  Signal  Processing,  Vol.  ASSP- 
32,  No.  2,  April,  1984. 

3.  Knapp,  C.  H.,  and  Carter  G.  C.,  “The  Generalized  Correlation  Method  for  Estimation 
of  Time  Delay”,  IEEE  Transactions  on  Acoustics,  Speech  and  Signal  Processing,  Vol. 
ASSP-24,  No.  4,  August,  1976. 

4.  Stein,  Seymour,  "Algorithms  for  Ambiguity  Function  Processing",  IEEE  Transactions 
on  Acoustics,  Speech  and  Signal  Processing,  Vol.  ASSP-29,  No.  3,  June  1981. 

5.  Brown,  W.  A.,  and  Spooner,  C.,  “Cyclostationary  Time  Difference  of  Arrival 
Techniques  for  Signals  Intelligence:  Part  I  of  H”,  Statistical  Signal  Processing  Inc.  and 
Mission  Research  Corporation  Inc.,  technical  report,  June  1996. 

6.  Gardner,  William  A.,  "Exploitation  of  Spectral  Redundancy  in  Cyclostationary 
Signals",  IEEE  Signal  Processing  Magazine,  Apn\,  1991. 

7.  Gardner,  William  A.,  ed.,  Cyclostationarity  in  Communications  and  Signal 
Processing,  IEEE  Press,  1994. 

8.  Chen,  Chih-Kang  and  Gardner,  William  A.,  "Signal-Selective  Time-Difference- 
of-  Arrival  Estimation  for  Passive  Location  of  Man-Made  Signal  Sourees  in 
Highly  Corruptive  Environments,  Part  I:  Theory  and  Method",  IEEE  Transactions 
on  Signal  Processing,  Vol.  40,  Number  5,  May  1992. 

9.  Chen,  Chih-Kang  and  Gardner,  William  A.,  "Signal-Selective  Time- 
Difference-of-  Arrival  Estimation  for  Passive  Location  of  Man-Made  Signal 
Sources  in  Highly  Corruptive  Environments,  Part  H:  Algorithms  and 
Performance",  IEEE  Transactions  on  Signal  Processing,  Vol.  40,  Number  5,  May 
1992. 

10.  Brown,  W.  A.,  Introduction  to  Second-Order  Cyclostationary  Signals,  Mission 
Research  Corp.  briefing,  April  22,  1996. 

11.  Schau,  H.  C.,  and  Robinson,  A.  Z.,  “Passive  Source  Location  Employing 
Spherical  Surfaces  from  Time-of-Arrival  Differenees”,  IEEE  Transactions  on 
Acoustics,  Speech  and  Signal  Processing,  Vol.  ASSP-35,  No.  8,  August  1987. 

12.  Loomis,  Herschel  H.  Jr.,  "Geolocation  of  Electromagnetic  Emitters", 

Unpublished  work.  Naval  Postgraduate  Sehool,  Monterey,  CA,  1993. 


97 


13.  Ho,  K.  C.  and  Chen,  Y.  T.,  “Solution  and  Performance  Analysis  of 
Geolocation  by  TDOA”,  IEEE  Transactions  on  Aerospace  and  Electronic 
Systems,  Vol.  29,  No.  4,  October  1993. 

14.  Rusu,  Petre,  "The  Equivalence  of  TOA  and  TDOA  RF  Transmitter  Location", 
Unpublished  work.  Applied  Research  Laboratory:  The  University  of  Texas  at 
Austin,  Austin,  TX,  1995. 

15.  Leach,  M.  and  Feuerstein,  D.,  “Preliminary  Test  Results  from  a  Carry-on 
Multi-platform  GPS-Assisted  Time  Difference  of  Arrival  System”,  Proceedings 
of  the  Symposium  on  Radiolocation  and  Direction  Finding,  San  Antonio,  Texas, 

June  1995. 

16.  Rusu,  P.  and  Giulianelli,  L.,  “A  Stochastic  Model  for  Time  Difference  of 
Arrival  Error  Estimation”, '^ARL:UT  Inter-office  Memorandum,  1995. 

17.  Applied  Research  Laboritories:  University  of  Texas  at  Austin,  briefing,  September 
1996. 

18.  Rusu,  Petre,  “On  the  Dilution  of  Precision  Factors  in  TDOA  Computations”, 
ARL:UT  Inter-office  Memorandum,  1995. 

19.  Loomis,  H.  H.  Jr.  And  Bernstein,  R.  F.  Jr.,  “Realization  of  TDOA  Estimation 
Architectures”,  IEEE  Transactions  on  Signal  Processing,  Vol  50,  No.  3,  August  1995. 

20.  Brown,  W.  A.  and  Loomis,  H.  H.  Jr.,  “Digital  Implementation  of  Spectral 
Correlation  Analyzers”,  IEEE  Transactions  on  Signal  Processing,  Vol.  41,  No.  2, 
February  1993. 

21.  Brown,  W.  A.  and  Loomis,  H.  H.  Jr.,  “A  Review  of  Digital  Spectral  Correlation 
Analysis:  Theory  and  Implementation”,  Cyclostationarity  in  Communications  and  Signal 
Processing,  William  A.  Gardner,  ed.,  IEEE  Press,  1994. 

22.  Roberts,  R.,  Brown,  W.  A.  and  Loomis,  H.  H.  Jr.,  “Computationally  Efficient 
Algorithms  for  Cyclic  Spectral  Analysis”,  IEEE  Signal  Processing  Magazine,  Vol  39 
No.  4,  April  1991. 

23.  MATIAB  Reference  Guide,  Version  4.2c,  The  Mathworks,  Inc.,  August  1992. 

24.  Simulink  User’s  Guide,  The  Mathworks,  Inc.,  1995. 


98 


BIBLIOGRAPHY 


Janiczek,  P.  M.,  ed.,  Global  Positioning  System  Volume  I,  The  Institute  of  Navigation, 
Washington,  D.C.,  1980. 

Leick,  Alfred,  GPS  Satellite  Surveying,  John-Wiley  &  Sons,  New  York,  NY,  1990. 

Parkinson,  Bradford  W.  and  Spilker,  James  J.  Jr.,  ed..  Global  Positioning  System: 

Theory  and  Applications  Volume  I,  American  Institute  of  Aeronautics  and  Astronautics, 
Inc.,  Washington,  D.C.,  1996. 

Parkinson,  Bradford  W.  and  Spilker,  James  J.  Jr.,  ed.,  Global  Positioning  System: 

Theory  and  Applications  Volume  II,  American  Institute  of  Aeronautics  and  Astronautics, 
Inc.,  Washington,  D.C.,  1996. 

Price,  Michael  G.,  "Mathematics  of  Geolocation",  Informal  report,  Systeka,  Inc., 
Seabrook,  MD,  1995. 

r 

Rusu,  Petre,  "The  Equivalence  of  TOA  and  TDOA  RF  Transmitter  Location", 
Unpublished  work.  Applied  Research  Laboratory:  The  University  of  Texas  at  Austin, 
Austin,  TX,  1995. 

Rusu,  Petre,  "A  TOA-FOA  Approach  to  Find  the  Status  of  RF  Transmitters", 
Unpublished  work.  Applied  Research  Laboratory:  The  University  of  Texas  at  Austin, 
Austin,  TX,  1995. 

Spooner,  Chad  M.,  "An  Overview  of  Recent  Developments  in  Cyclostationary  Signal 
Processing",  Proceedings  of  the  Second  Workshop  on  Cyclostationary  Signals,  Mission 
Research  Corporation,  Monterey,  CA,  August  1-2,  1994. 


99 


100 


APPENDIX  A.  MATLAB  CODE  FOR  ALGORITHMS 


The  programs  that  follow  constitute  the  major  sections  of  code  used  in  this  work. 
They  all  contain  headers  and  commenting  for  explanation  purposes.  They  are  described 
briefly  in  Table  A-1  below. 


Program  Title 

Purpose 

loaddata.m 

Loads  headers  and  signal  data,  converts  4  bytes  real,  4 

bytes  imaginary  to  complex  vector  notation. 

ssca.m 

Strip-spectral  correlation  analyzer.  Calculates  the  SCF 

of  the  input  sequence  along  diagonal  lines  according  to 

the  SSCA  algorithm. 

plotscf.m 

Transforms  the  diagonal  lines  of  ssca.m  into  lines  along 

f  and  a.  Plots  both  the  SCF  and  the  magnitude  of  SCF 

vs  cycle  frequency. 

speccoa.m 

Resamples  time  domain  data  to  user  specified  value  (p). 

Calculates  the  TDOA  of  two  input  sequences  in  terms 

of  samples  using  SPECCOA  algorithm. 

rusugeo.m 

Takes  positions  of  observers  and  three  TDOA(s)  to 

produce  geolocation  in  WGS84  ECEF  coordinates. 

Finds  both  solutions.  User  must  eliminate  one. 

Table  A-1  -  Summary  of  code 


^Hi********Hc******Hi******************Hi*****5k************************ 

% 

%  Name:  David  A  Streight 
% 

%  Naval  Postgraduate  School  -  Monterey,  California 
% 

%  loaddata.m 
% 

%  Operating  Environment 
%  Windows  95 
% 

%  Description 

%  Loads  32-bit  float  ARL:UT  data  into  four  header  and  data  vectors. 

%  Master  site  is  obsl,  other  sites  are  2,  3  and  4. 

% 

%  Date  of  last  revision 
%  20  December  1996  ' 

% 

%  Inputs 

%  Four  ARL.'UT  stored  files 

% 

%  Outputs 
%  None 

% 

<^**!|!***  *****!|:**!C**H:********N:*!t:*!)!!C**!l:!|:!|:*********H!  iic^H:*!)!****!!:**  ******* 
^**^:******H:********************H:**H:*******H:***:|:*  ********  **********:!: 

clear  all; 

%  open  four  observer  files  read  only 

obsl  =  fopen('F:\data\nov_ut~2\mcc\409400.revVr'); 
obs2  =  fopen('F:\data\nov_ut~2\berg\409400.revVr'); 
obs3  =  fopen('F:\data\nov_ut~2\beecaves\409400.revVr'); 
obs4  =  fopen('F:\data\nov_ut~2\expo\409400.rev','r’); 

%  load  3  header  variables  into  header  vectors 
[header  1]  =  fread(obs  1,3, 'float'); 

[header2]  =  fread(obs2,3, 'float'); 

[header3]  =  fread(obs3,3,'float'); 

[header4]  =  fread(obs4,3, 'float'); 

%  load  data  into  4  data  matrices  col#l  -  real  col#2  -  imaginary 
[xl]  =  fread(obsl,[8192,2],’float'); 

[x2]  =  fread(obs2,[8192,2],'float'); 

[x3]  =  fread(obs3,[8192,2],'float'); 

[x4]  =  fread(obs4,[8192,2],'float’); 
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%  convert  data  matrices  to  complex  vectors 
si  =xl(:,l)+j*xl(:,2); 
s2  =  x2(:,l)+j*x2(:,2); 
s3  =  x3(:,l)+j*x3(:,2); 
s4  =  x4(:,l)+j*x4(:,2); 

%  clear  xi’s 
clear  xl; 
clear  x2; 
clear  x3; 
clear  x4; 

%  close  all  files 
fclose('air); 


103 


^***H:****!|:*!)!*=|:*!!=**>t:!|:!t:*H:H:***!t:***si!*********!|:*!)!%*H:*!|!**!i!!(:***5l:H:**j(c5k4:s|:*** 

% 

%  Name:  David  A  Streight 
% 

%  Naval  Postgraduate  School  -  Monterey,  California 
% 

%  ssca.m 
% 

%  Operating  Environment 
%  Windows  95 
% 

%  Description 

%  Calculates  the  auto  spectral  density  function  using  the  strip 
%  spectral  correlation  algorithm  (ssca) 

% 

%  Date  of  last  revision 
%  18  November  1996 
% 

%  Inputs 
%  x-SOI 

% 

%  Outputs 

%  SDF  -  the  results  of  the  ssca  computation.  Must  be  converted 
%  and  plotted  by  plotsdf.m 
% 


x  =  [sl(l;2^12,l)].’; 

%  prepare  variables 

s  =  size(xj; 

N  =  s(1,2); 

%  number  of  points  in  vector 

Np  =  8; 

%  course  FFT  number 

xp  =  [x  zeros(l,Np)]; 

%  zero  pad  input  variable  for  Np  length  slices 

w  =  (hamming(Np))'; 

%  hamming  window  for  course  FFT 

W  =  (haniming(N))'; 

%  hamming  window  for  P  point  FFT  at  the  end 

L  =  Np/4; 

%  decimation  factor 

P  =  floor(N/L); 

p  =  zeros(l,floor(P/L)); 

xT  =  zeros(l,Np); 

k  =  -Np/2:Np/2-l; 

dm  =  zeros(l,Np); 

XTi  =  zeros(N,Np); 

SXT  =  zeros(Np,N); 
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%  find  complex  demodulates 
forp  =  0:P-l 

xT  =  fft((xp(l,p*L+l:p*L+Np)  .*  w),  Np);  %  take  Np  point  FFT  of  windowed  data 
xT  =  [xT(l,5;8)  xT(l,l:4)];  %  fold  fft  to  account  for  neg  freq's 

dm  =  exp(-j*2*pi*k*p*L/Np);  %  calculate  freq  shifts  for  this  set 

XTi(p*L+l ,:)  =  xT  .*  dm;  %  freq  shift  and  populate  odd  col's  of  interp 

matrix 
end 

%  interpolate  ***works  for  Np  =  8  only  at  this  point*** 
for  V  =  1:P-1 

XTi(L*v,:)  =  (XTi(L*v- 1 ,;)  +  XTi(L*v+l  ,:))./2;  %  fill  in  even  col’s  linearly 

end 

XTi(N,:)  =  XTi(  1 %  hold  next  to  last  col  for  last  col  values 
%  finish  up 

xm  =  [x;x;x;x;x;x;x;x]';  %  matrix  of  original  signal 

% Wm  =  [W ;W ;W ; W ;W ;W ; W ;W] ;  %  matrix  of  column-wise  hamming 

window 

SXT  =  fft(XTi  .*  xm,  N);  %  take  N-point  FFT  of  products 

SDF  =  SXT.';  %  transpose  after  FFT 

SDF  =  [SDF(:,1  :P- 1)  SDF(:,P:N)];  %  fold  fft  to  account  for  negative  freq's 
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<%,******:|=*****=|:*=)=*****=|r*=k*******H=***H=**=l=**********=C*******H:H=H=H:****** 

%********H=*Hc******H=**************=^*H=*****************:.=  *****:^******* 
% 

%  Name:  David  A  Straight 
% 

%  Naval  Postgraduate  School  -  Monterey,  California 
% 

%  PLOTSDF:  Convert  and  plot  SSCA  output 
% 

%  Operating  Environment 
%  Windows  95 
% 

%  Description 
% 

% 

% 

% 

%  Date  of  last  revision 
%  15  November  1996 
% 

%  Inputs 

%  Matrix  of  Spectral  Density  Function  values  from  the  SSCA  algorithm 
%  SDF 

% 

%  All  variables  used  in  the  SSCA  computation  (no  clear  executed) 

% 

%  Outputs 

%  Plot  in  3D  of  the  SDF  function  drawn  along  lines  of  cycle  frequency 
^^****=l=*Hc=t=***=H*******H=*********=l=**********=)=********=l=  =(=****=(:****:(=***** 

%******:,:***********=(=  H=*****Hc*H=***=C**=i=**=l=*=t=:f=*:^=i=*H=****H=:f:**%X:******4=*=tc* 


.%  determine  the  span  of  alpha 
alphamax  =  2*N  -  2*N/Np  -  2; 
alphamin  =  -2*N; 
alphaspan  =  alphamax  -  alphamin; 

%  setup  intermediate  matrices 
I  =  zeros(alphaspan+l,Np); 

X  =  zeros(alphaspan+l,Np); 

Y  =  zeros(alphaspan+l,Np); 


% 

maxalpha  =  zeros(l,alphamax+l); 

%  populate  matrices  from  SDF 
for  q  =  -N/2:N/2-l 
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for  k  =  -Np/2:Np/2  -  1 

I(2*N*k/Np+2*q+2*N+l,k+Np/2+l)  =  abs(SDF(k+Np/2+l,q+N/2+l)); 
end 

end 

%  normalize  magnitudes 
I  =  I/max(max(I)); 

%  populate  f  and  alpha  coordinate  matrices 

j=l:8; 

y  =  ones(l,8); 

for  i  =  l:alphaspan+l 

X(i,:)  =  2*N*G-Np/2-l)/Np-(i-2*N-l)/2; 

Y(i,:)  =  ((i-2*N-l).*y)/2; 


end 

%  plot  the  positive  alpha  half  of  the  plane 
figure(l); 

surf(X(alphaspan-alphamax :  alphaspan+ 1 , ;  ),Y(alphaspan- 
alphamax:alphaspan+l  ,:),I(alphaspan-alphamax:alphaspan+l 
view(-105,60),grid,axis([-20000  10000  0  8000  0  1]) 

%  now  plot  alpha  vs  mag  of  SDF  for  those  values 
for  1  =  l;alphamax+l 

maxa(l)  =  max(I(alphaspan-alphamax+l-l,:)); 
end 

maxalphas  =  nonzeros(maxa); 
maxalphas  =  maxalphas'; 
sz  =  size(maxalphas); 
alphanum  =  0:sz(l,2)-l; 
figure(2); 

plot(alphanum,maxalphas) 

clear  i; 
clear  j; 
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^:i:5f::t:******************  ********************************  ****:!:******** 

<^HJ****  **************************  Hi***  ******************  :i;*5f:*  *:}:******* 
% 

%  Name:  David  A  Streight 
% 

%  Naval  Postgraduate  School  -  Monterey,  California 
% 

%  SPECCOA  -  Cyclostationary  Time  Difference  of  Arrival  Algorithm 
% 

%  Operating  Environment 
%  Windows  95 
% 

%  Description 

%  Produces  an  estimate  of  the  time  difference  of  arrival  of  a 
%  signal  of  interest  between  two  spatially  separated  receivers 
%  in  terms  of  multiples  of  the  sampling  frequency 
% 

%  Date  of  last  revision 
%  06  January  1997 

% 

%  Inputs 

%  p  -  up  sample  integer 
%  q  -  down  sample  integer 

%  terms  -  number  of  terms  used  in  linear  interpolation 
%  N  -  number  of  samples  of  decimated  sequence  to  use 
%  a  -  cycle  frequency  of  interest  in  terms  of  decimated  fs/N 
%  soil  -  signal  captured  at  receiver  1  (N  samples) 

%  SOI2  -  signal  captured  at  receiver  2  (N  samples) 

%  *both  SOI(s)  sampled  at  rate  fs  samples  per  sec 
% 

%  Outputs 

%  Plot  of  SPECCOA  function  versus  samples.  Peak  is  TDOA  estimate 
% 

^******************************:j::{::j;:ij:}::^:{::fi:{::{::j::i::|::|::jc:^:j.:j;:j.:j:^^^^:j;:^^^^:^^5jj^^^jj, 

^*******************************:^c:f::i::{::}c:};:i;:}::ic:{::i::f;:{::{:5j::j;^:j;:j:^;j:;j;^:|;^^jj;^:|;^^5j;^^^ 


loaddata; 

%  load  decimated  signals 

%  resample  variables 

p  =  256; 

%  up  sample 

q=l; 

%  down  sample 

terms  =  3; 

%  terms  in  the  linear  interpolation 

I  =  p/q; 

%  total  interpolation 

%  cyclostationary  variables 

a  =  2049  *  I; 

%  cycle  frequency  at  new  sample  rate 
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N  =  4096  *  I; 
fs  =  9765.625; 

%  prepare  loop  variables 
SXT  =  zeros(l,N); 
SYXT  =  zeros(l,N); 


%  length  of  sequence  at  new  sample  rate 


%SCFofSOIl  atalphaO 
%  Cross  SCF  at  alphaO 


Sigl=sl(l:(N/I),l).'; 

Sig2  =  s2(l:(N/I),l).’; 

soil  =  resample(Sigl,p,q, terms); 
SOI2  =  resample(Sig2,p,q,terms); 

Xf  =  fft(SOIl,N); 

XT  =  fftshift(Xf); 

Yf  =  fft(SOI2,N); 

YT  =  fftshift(Yf); 


%  set  up  vectors 


%  resample  and  set  as  SOI! 
%  resample  and  set  as  SOI2 

%  freq  spectrum  of  SOU 
%  adjust  for  MATLAB  fft 
%  freq  spectrum  of  SOI2 
%  adjust  as  above 


XT  =  [zeros(l,a/2)  XT  zeros(l,a/2)];  %  zero  pad  at  each  end  for  correlation  loop 
YT  =  [zeros(l,a/2)  YT  zeros (1, a/2)];  %  ditto 


SXT(1,:)  =  XT(1,1:N)  .*  conj(XT(l,a+l:N+a)); 
SYXT(1,:)  =  YT(1,1:N)  .*  conj(XT(l,a+l:N+a)); 


isp  =  ifft((conj(SXT)  .*  SYXT),  N); 

s  =  -N/2;N/2-l; 

kernel  =  exp(j*pi*a*s/N); 

sp  =  fftshift(isp)  .*  kernel; 

plot(s,real(sp)) 


%  IFFT  of  complex  product  from  above 
%  time  in  samples 
%  freq  kernel  at  alphaO 
%  adjust  and  multiply 
%  plot  to  find  TDOA 
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^  :ic  H:  ^  ^ ^  Hi  if:  4c  sK  ^  ^  Hi  ^  ^  Hi  5k  %  ^  ^  Hi  ^  Hi  Hi  ^  =t:  =1:  ^  :f:  sj:  sN  ^  ij: 5i;  % 

Hi  Hi  *  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi*  Hi  Hi  ♦  Hi  Hi  Hi  Hi  **  Hi  Hi  Hi  Hi  Hi  Hi  Hi  *  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  Hi  *  Hi  Hi  Hi  Hi  Hi  Hi  Hi 

% 

%  Name:  David  A  Streight 
% 

%  Naval  Postgraduate  School  -  Monterey,  California 
% 

%  Closed  Form  TDOA  Geolocation  Solution 
% 

%  Operating  Environment 
%  Windows  95 
% 

%  Description 

%  Produces  two  WGS84  based  3-dimensional  geolocation  solutions  for  a  transmitter 
%  given  the  location  of  each  of  four  observers  and  the  TDOAs  for  their  observations 
%  Data  outside  the  algorithm  must  be  used  to  eliminate  the  ambiguity. 

% 

%  Date  of  last  revision 
%  10  October  1996 

% 

%  Inputs 

%  WGS84  based  x,  y  and  z  coordinates  for  each  of  four  observers 
%  xi,  yi  zi,  i  =  1,2,3,4 
% 

%  TDOAs  for  the  observers 

%  Dij,  i  =  1, 2,3,4  and  j  =  1, 2,3,4  where  i  *  j  and  Dij  =  -Dji 

% 

%  Outputs 

%  WGS84  based  x,  y  and  z  coordinates  for  the  transmitter 
%  rtl  and  rt2  (vectors) 

% 

^*:i:************  ************************************  **5ii** 

^  *******************:}::*:  :jc :}:  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  jjj  ^  ^  jjj  jj.  jj.  jj.  jj,  jj,  ^  ^  jj. 


%  positions 

xt  =  -741941.47;  yt  =  -5462120.41;  zt  =  3198242.728; 
xl  =  -741203.52;  yl  =  -5456323.70;  zl  =  3208396.55; 
x2  =  -735740.31;  y2  =  -5467456.78;  z2  =  3190514.57; 
x3  =  -754177.92;  y3  =  -5459066.31;  z3  =  3200764.99; 
x4  =  -731248.92;  y4  =  -5463237.87;  z4  =  3198800.64; 

%  set  up  inputs  and  utility  vectors 

N  =  0;  %  relative  times  of  arrival  using  observer  1  as 

baseline 

t2  =  -D12;  %  D12  is  defined  as  tl  - 12  thus  t2  =  -D12  if  tl  =  0 

t3  =  -D13;  %  ditto 
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%  ditto 


t4  =  -D14; 

rl  =  [xl;yl;zl];  %  observer  position  vectors 

t2  =  [x2;y2;z2]; 
r3  =  [x3;y3;z3]; 
r4  =  [x4;y4;z4]; 

Euclid_sq  1  =  x  1  -^2  +  y  1  ^^2  +  z  1  ^^2;  %  squares  of  Euclidean  distances 

Euclid_sq2  =  x2'^2  +  y2^2  +  z2^2; 

Euclid_sq3  =  x3^2  +  y3^2  +  z3^2; 

Euclid_sq4  =  x4^2  +  y4^2  +  z4'^2; 

%  set  up  solution  matrices 

A  =  [xl-x2  yl-y2  zl-z2;x2-x3  y2-y3  z2-z3;x3-x4  y3-y4  z3-z4]; 

q  =  [tl-t2;t2-t3;t3-t4]; 

s  =  zeros(3,l);  *' 

s(l,l)  =  (Euclid_sql  -  Euclid_sq2  -  tl'^2  + 12^2).*0.5; 
s(2,l)  =  (Euclid_sq2  -  Euclid_sq3  - 12^2  +  t3''2).*0.5; 
s(3,l)  =  (Euclid_sq3  -  Euclid_sq4  - 13^2  + 14^2).*0.5; 

g  =  A\q;  %  use  matrix  division  (MATLAB)... 

%  in  lieu  of  the  inverse  function... 
h  =  A\s;  %  to  solve  this  part 

%  find  terms  for  at^2  +  bt  +  c 

a  =  ((g(El))'^2  +  (g(2,l))''2  +  (g(3,l))^2)  -  1;%  find  quadratic  term 
b  =  (g.'  *  h)  +  (g.'  *  rl)  +  tl ;  %  find  first  order  term 

c  =  (norm(h  -  rl).^2)  -  tl^2;  %  find  constant 

p  =  [a  2*b  c];  %  form  polynomial  in  t 

%  find  the  two  times  as  roots  of  this  equation 

t  =  roots(p);  %  find  roots  of  at‘^2  +  2bt  +  c 

%  substitute  into  t-parametrized  equation  for  r  and  solve  for  real  part  of  two  solutions 
rtl  =real(g.*t(l,l)  +  h); 
rt2  =  real(g  .*  t(2,l)  +  h); 
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APPENDIX  B.  PLOTS  FOR  ARL:UT  TEST  DATA 


The  following  150  figures  are  the  result  of  testing  25  the  cyclostationary  TDOA 
processor  for  25  time  periods.  For  each  time  period,  6  figures  appear:  one  SCF  plot  for 
each  of  the  four  observers  grouped  two  to  a  figure,  one  magnitude  of  SCF  vs  cycle 
frequency  for  each  of  the  four  observers  grouped  two  to  a  figure  and  three  SPECCOA 
plots  grouped  two  and  one  to  a  figure. 

For  the  SCF  plots,  only  the  positive  cycle  frequency  values  are  plotted.  The 
spectral  frequency  axis  runs  from  the  upper  left  to  the  lower  right;  the  cycle  frequency 
axis  runs  from  the  upper  right  to  the  lower  left.  Both  the  spectral  frequency  and  the 
cycle  frequency  are  plotted  in  terms  of  fj/N  where  fj  =  9765.625  Hz  and  N  =  4096.  The 
magnitude  is  of  course  the  axis  pointing  to  the  top  of  the  page.  The  magnitude  of  the 
plots  is  nomalized  by  the  maximum  value  of  the  SCF  and  thus  ranges  from  0  to  1. 

The  magnitude  of  the  SCF  versus  cycle  frequency  plots  also  display  only  the 
positive  cycle  frequencies.  The  horizontal  axis  is  cycle  frequency  in  terms  of  fj/N  as 
above.  The  vertical  axis  is  again  normalized  magnitude.  These  plots  merely  represent 
looking  at  the  SCF  plots  from  the  cycle  frequency  axis  side.  They  aided  in  choosing  the 
best  cycle  frequency  for  SPECCOA. 

Finally,  the  SPECCOA  plots  show  the  output  of  the  speccoa  code.  The 
horizontal  axis  represents  the  TDOA  in  terms  of  the  numbers  of  samples  of  delay.  The 
horizontal  axis  is  once  again  magnitude. 
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Figure  B-6  -  SPECCOA  for  master  and  slave  #3,  time  409400 


Figure  B-10  -  Cycle  freq  vs  Mj 


Figure  B-12  -  SPECCOA  for  master  and  slave  #3,  time  409405 
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Figure  B-15  -  Cycle  freq  vs  Max-magnitude  of  SCF  for  master  (top)  and  slave  #1  (bottom),  time  409410 
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Figure  B-18  -  SPECCOA  for  master  and  slave  #3,  time  409410 
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Figure  B-21  -  Cycle  freq  vs  Max-magnitude  of  SCF  for  master  (top)  and  slave  #1  (bottom),  time  409415 
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Figure  B-24  -  SPECCOA  for  master  and  slave  #3,  time  409415 
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Figure  B-27  -  Cycle  fn 


Figure  B-28  -  Cycle 


Figure  B-33  -  Cycle  fireq  vs  Max-magnitude  of  SCF  for  master  (top)  and  slave  #1  (bottom),  time  4094 
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Figure  B-34 


-  Cycle  freq  vs  Max-magnitude  of  SCF  for  slave  #2  (top)  and  slave  #3  (bottom), 


time 
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Figure  B-36  -  SPECCOA  for  master  and  slave  #3,  time  409425 
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Figure  B-39  -  Cycle  fr( 
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Figure  B-54  -  SPECCOA  for  master  and  slave  #3,  time  409455 
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Figure  B-58  >  Cycle  freq  vs  Mj 
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Figure  B-60  -  SPECCOA  for  master  and  slave  #3,  time  409460 
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Figure  B-63  -  Cycle  freq  vs  Max-magnitude  of  SCF  for  master  (top)  and  slave  #1  (bottom),  time  409465 
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Figure  B-66  -  SPECCOA  for  master  and  slave  #3,  time  409465 
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Figure  B-70  -  Cycle  freq  vs  M: 
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Figure  B-78  -  SPECCOA  for  master  and  slave  #3,  time  409475 
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Figure  B-81  -  Cycle  freq 
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Figure  B-84  -  SPECCOA  for  master  and  slave  #3,  time  409480 


Figure  B-87  -  Cycle  freq  vs  Max-magnitude  of  SCF 


Figure  B-90  -  SPECCOA  for  master  and  slave  #3,  time  409485 
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Figure  B-96  -  SPECCOA  for  master  and  slave  #3,  time  409490 
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Figure  B-99  -  Cycle  fireq  vs  Max-magnitude  of  SCF  for  master  (top)  and  slave  #1  (bottom),  time  409495 
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Figure  B-102  -  SPECCOA  for  master  and  slave  #3,  time  409495 
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Figure  B-105  -  Cycle  fireq  vs  Max-magnitude  of  SCF  for  master  (top)  and  slave  #1  (bottom),  time  409500 
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Figure  B-106  -  Cycle  freq  vs  Max-magnitude  of  SCF 
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Figure  B-108  -  SPECCOA  for  master  and  slave  #3,  time  409505 
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Figure  B-1 1 1  -  Cycle  freq  vs  Max-magnitude  of  SCF  for  master  (top)  and  slave  #1  (bottom),  tim« 
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Figure  B-1 14  -  SPECCOA  for  master  and  slave  #3,  time  409505 
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Figure  B-1 17  -  Cycle  freq  vs  Max-magnitude  of  SCF  for  master  (top)  and  slave  #1  (bottom),  ti; 
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Figure  B-1 18  -  Cycle  freq  vs  Max-magnitude  of  SCF  for  slave  #2  (top)  and  slave  #3  (bottom),  tin 
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Figure  B-126  -  SPECCOA  for  master  and  slave  #3,  time  409515 


Figure  B-129  -  Cycle  freq  vs  Max-magnitude  of  SCF  fc 
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Figure  B-130  -  Cycle  freq  vs  Max-magnitude  of  SCF  for  slave  #2  (top)  and  slave  #3  (bottom),  time  409520 
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Figure  B-132  -  SPECCOA  for  master  and  slave  #3,  time  409520 
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Figure  B-136  -  Cycle  freq  vs  Max-magnitude  of  SCF  for  slave  #2  (top)  and  slave  #3  (bottom),  ti 
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Figure  B-138  -  SPECCOA  for  master  and  slave  #3,  time  409525 
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Figure  B-141  -  Cycle  freq  vs  Max-magnitude  of  SCF  for  master  (top)  and  slave  #1  (bottom),  time  409530 
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Figure  B-144  -  SPECCOA  for  master  and  slave  #3,  time  409530 
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Figure  B-147  -  Cycle  freq  vs  Max-magnitude  of  SCF  for  master  (top)  and  slave  #1  (botto 
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Figure  B-150  -  SPECCOA  for  master  and  slave  #3,  time  409535 
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